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Abstract 

To  aid  the  viewer  in  interpreting  the  information  contained  in  an  image,  spati2d  resolution  of 
the  image  can  be  increased.  One  way  to  increase  the  resolution  is  to  improve  the  optics  and  sensors 
of  the  imaging  system.  Another  is  to  process  the  data  on  the  ground  by  computer.  Enhancement 
techniques  caa  be  employed  to  alter  the  image  and  extract  necessary  data.  Current  techniques  to 
manipulate  and  enlarge  an  image  are  nearest-neighbor,  bilinear  interpolation,  and  cubic  convolution 
which  can  give  greater  clarity  and  readability.  To  satisfy  the  growing  demand  for  quality  (high- 
resolution)  remotely  sensed  images,  new  techniques  to  enhance  images  need  to  be  developed. 

Kriging  is  the  geostatistical  spatial  estimation  technique  developed  by  Matheron  that  opti¬ 
mally  predicts  spatial  data  using  observations  at  nearby  locations.  Kriging  can  be  applied  to  many 
specialties  and  applications  such  as  lead  contamination  in  the  soil  near  a  smelter  or  prevalence  of 
AIDS  within  a  population  viewed  geographically. 

This  research  investigates  the  application  of  kriging  as  a  technique  of  digital  image  resolu¬ 
tion  enhancement  by  adapting  its  use  for  the  personnel  computer  and  validating  its  capability 
and  reliability.  Twenty-seven  percent  of  the  images  tested  produced  significantly  lower  means  for 
kriging  than  cubic  convolution  with  45  percent  of  the  kriged  images  producing  lower  variances. 
When  universal  and  ordinary  kriging  were  tested,  no  difference  was  found  between  the  means  and 
variances  while  no  significant  differences  between  the  use  of  a  model  and  lineeir  interpolation  for 
semi-variogram  estimation  was  found. 
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ANALYSIS  OF  KRIGING  APPLIED  TO 
RESOLUTION  ENHANCEMENT  OF 
DIGITAL  SATELLITE  IMAGERY 

I.  Introduction 

This  research  continues  the  efforts  of  Captain  Donald  W.  McGee,  Jr.  and  investigates  the 
application  of  statistical  kriging  to  digital  image  resolution  enhancement.  This  research  builds 
on  the  theoretical  background  of  McGee  and  further  adapts  the  technique  of  kriging  by  taking 
advantage  of  the  regular  nature  of  satellite  imagery  to  simplify  processing. 

Kriging  is  the  geostatistical  spatial  estimation  technique  developed  by  Matheron  that  opti¬ 
mally  predicts  spatial  data  using  observations  at  nearby  locations  (18:1259).  Christensen  states 
“the  object  of  kriging  is  the  prediction  of  unobserved  spatial  random  variables  based  on  the  values 
of  observed  random  variables”  (4:267)  and  adds, “there  are  innumerable  situations  in  which  data 
are  collected  at  various  locations  in  space  and  thus  innumerable  potential  application  for  methods 
of  analysis  for  spatial  data”  (4:262).  Kriging  can  be  applied  to  many  specialties  and  applications 
such  as  lead  contamination  in  the  soil  near  a  smelter  or  prevalence  of  AIDS  within  a  population 
viewed  geographically  (4:263)  or  resolution  enhancement  of  digital  imagery  .  This  paper  provides 
background,  objectives,  summary  of  current  knowledge,  scope,  and  methodology  of  the  research  to 
apply  kriging  to  digital  image  resolution  enhancement. 


1.1  Background 

The  traditional  method  for  extracting  data  from  images  is  manual  interpretation  (3:719). 
With  the  invention  of  the  charge-coupled  device  (CCD)  and  advances  in  computer  technology, 
quick  and  easy  spatial  resolution  enhancement  is  now  possible.  Spatial  resolution  is  the  degree 
of  discernible  detail  of  an  image  or  the  ability  to  distinguish  between  spatial  elements  within  an 
image  (17:36).  Because  increasing  the  spatial  resolution  of  imagery  extracts  the  most  information 
possible  from  the  image,  it  has  become  major  concern  of  remote  sensing  specialists,  astronomers, 
and  the  intelligence  community.  Resolution  may  be  increased  in  two  ways: 

1.  Optics  and  sensors  of  the  imaging  system  can  be  improved  (larger  optics,  lower  altitudes, 
or  larger  CCD  arrays)  to  increase  the  spatial  resolution  of  the  images  produced.  Changing 


1-1 


the  imaging  system  is  undesirable  because  it  is  costly  and  produces  massive  amounts  of  data 
creating  problems  with  data  handling  and  distribution. 

2.  Image  data  can  be  processed  on  the  ground  to  enhance  resolution.  Processing  involves  enlarg¬ 
ing  the  original  image  (through  the  addition  of  blank  pixels  to  produce  a  larger  image)  and 
then  estimating  and  assigning  new  values  to  the  coordinates  of  the  blank  pixels  to  reproduce 
the  image.  This  enhancement  method  is  preferred  due  to  its  lower  cost  and  relative  simplicity. 

Though  processing  of  the  images  to  enhance  their  resolution  is  preferred,  problems  exist  with 
the  current  techniques.  Three  processing  methods  are  commonly  used  to  estimate  pixel  values;  they 
are  nearest-neighbor,  bilinear  interpolation,  and  cubic  convolution.  The  nearest-neighbor  technique 
leaves  images  blocky  in  appearance  while  bilinear  interpolation  tends  to  blur  the  image.  Neither 
result  is  desirable. 

This  research  will  continue  the  development  of  kriging  as  a  viable  and  preferred  image  res¬ 
olution  enhancement  technique.  On  kriging ’s  ability  of  spatial  estimation,  Cressie  states  “In  all 
comparisons,  on  real  and  simulated  data,  universal  kriging  generally  did  as  well  or  better  than  the 
other  methods”  (6:202).  Kriging  can  perform  as  well  or  better  than  current  enhancement  tech¬ 
niques  because  it  is  an  exact  method  of  estimation.  In  digital  image  terms,  the  gray  values  of 
the  original  image  are  reproduced  exactly  with  no  error  while  estimating  new  spatial  data.  Kriged 
images  will  not  have  the  blocky  appearance  of  nearest-neighbor  nor  blurs  the  image  as  does  bilinear 
interpolation  (20:112).  It  is  the  best  of  all  processing  techniques  and  is  truly  a  best  linear  unbiased 
predictor  (4:262). 


1.2  Research  Objective 

This  research  investigates  the  application  of  kriging  as  a  technique  of  digital  image  resolution 
enhancement  by  adapting  its  use  for  the  personnel  computer  and  validating  its  capability  and 
reliability. 

Ultimate  goals  of  the  research  are  streamlining  the  kriging  approach  (to  keep  processing  times 
low)  and  validation  of  previous  work  by  testing  the  technique  over  image  varying  in  terrain  and 
spectral  qualities. 


1.3  Summary  of  Current  Knowledge 

Many  techniques  of  generating  new  pixel  values  from  old  are  available  but  processing  costs 
and  quality  control  the  choice  of  method. 
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The  nearest-neighbor  technique  of  resolution  enhancement  is  the  simplest  method  in  use 
today.  The  value  of  the  new  pixel  takes  the  value  of  the  original  pixel  nearest  to  its  location.  The 
nearest-neighbor  technique  has  two  main  advantages;  it  is  fast  (takes  little  processor  time)  and 
preserves  the  exact  pixel  values  in  the  input  image  (17:138).  This  method  introduces  spatial  errors 
into  the  image  and  shifts  the  original  pixel  from  its  current  location  (3:734).  The  resultant  image 
has  a  blocky  appearance  and  can  be  difficult  to  read. 

Bilinear  interpolation  uses  the  gray  levels  of  the  four  nearest  neighbors  to  estimate  the  new 
pixel  values.  This  is  done  by  taking  the  gray  levels  of  the  four  nearest  neighbors  and  interpolating 
the  new  values  to  estimate  the  new  points.  This  technique  fills  the  desired  pixels  and  produces 
smoother  images  (13:301). 

Cubic  convolution  is  “based  on  the  fitting  of  a  two-dimensional,  third-degree  polynomial 
surface  to  the  region  surrounding  the  point”  (17:141)  of  the  form  (sina:)/r  and  produces  a  smooth 
estimate  of  the  gray  level  (13:300).  Mather  states,  this  technique  “gives  a  more  natural-looking 
image  without  the  blockiness  of  the  nearest- neighbor  [technique]  or  the  over-smoothing  of  the 
bilinear  method  [but  loses  some  data  in  the  process]”  (17:141).  One  problem  in  the  implementation 
of  this  method  is  that  it  requires  an  infinite  amount  of  data  to  accurately  estimate  the  (sinx)/x 
function  (2:55).  This  method  takes  more  processing  time  to  achieve  the  desired  results  than  the 
previous  methods  and  is  currently  the  “best”  technique  in  use  today. 

The  application  of  kriging  to  the  problem  of  digital  image  resolution  enhancement  has  produce 
promising  results.  McGee  states  “quantitative  comparison  of  the  kriging  and  cubic  convolution  tech¬ 
niques  indicates  that  image  enhancement  by  universal  punctual  kriging  can  produce  images  which 
are  more  faithful  than  cubic  convolution  to  the  original  image”  (19:4-30).  Kriging  produced  image 
differences  that  were  6.4  times  smaller  than  the  mean  difference  produced  by  cubic  convolution 
with  a  standard  deviation  of  4.73  times  smaller  (in  the  cases  examined)  (19).  “Kriging  methods 
were  clearly  superior  to  the  convolution  methods  in  recreating  the  original  image”  (19). 

Other  potential  applications  of  kriging  to  image  enhancement  suggested  by  McGee  include 
noise  reduction  and  as  a  method  of  edge  detection.  “Noise  elimination  using  kriging  should  produce 
better  results  than  other  methods  and  should  improve  image  classification  and  the  forecasts  made 
for  the  enhanced  imagery”  (19:5-2).  When  compared  to  the  original  image,  the  concentration  of 
error  produced  at  the  perimeter  of  objects  within  an  image  suggests  that  kriging  has  potential  as 
a  method  of  edge  detection  for  pattern  recognition  (19:5-2). 
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1-4  Scope 


This  research  will  be  limited  to  the  application  of  kriging  to  the  enhancement  of  image 
resolution.  Current  techniques  to  enhance  image  resolution  such  as  nearest-neighbor,  bilinear  in¬ 
terpolation,  and  cubic  convolution  will  be  compared  to  kriging.  Validation,  consolidation,  and 
streamlining  of  processing  will  be  the  main  emphasis  of  the  approach.  The  statistical  background 
of  kriging  and  application  of  kriging  for  data  reduction,  noise  reduction,  or  edge  detection  will  not 
be  covered  within  the  scope  of  this  analysis.  Due  to  technique  limitations,  visual  image  will  not  be 
included  in  the  research  with  the  statistical  analysis  be  the  limit  of  the  results. 
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II.  Literature  Review 


2. 1  Introduction 

Remote  sensing  is  the  measurement  and  collection  of  spectral  data  (reflected  or  emitted 
electromagnetic  (EM)  radiation)  from  a  planet’s  surface  or  atmosphere  from  a  remote  distance 
(17;1).  The  images  that  are  produced  from  this  data  are  representations  that  vary  in  brightness, 
color,  or  reflectance  and  can  be  described  mathematically  by  a  function  of  two  variables  (14:9). 
This  two-dimensional  light-intensity  function  f{x,y)  gives  the  brightness  of  the  image  (13:28). 

To  aid  the  viewer  in  interpreting  the  information  contained  in  an  image,  spatial  resolution  of 
the  image  can  be  increaised.  One  way  to  increase  the  resolution  is  to  improve  the  optics  and  sensors 
of  the  imaging  system.  Another  is  to  process  the  data  on  the  ground  by  computer.  Enhancement 
techniques  can  be  employed  to  alter  the  image  and  extract  necessary  data.  Current  techniques  to 
manipulate  and  enlarge  an  image  are  nearest-neighbor,  bilinear  interpolation,  and  cubic  convolution 
which  can  give  greater  clarity  and  readability.  These  techniques  have  their  problems  though,  either 
leaving  the  image  blocky  in  appearance  or  blur  the  image  (20:112). 

To  satisfy  the  growing  demand  for  quality  (high-resolution)  remotely  sensed  images,  new 
techniques  need  to  be  developed.  According  to  McGee,  the  geostatistical  estimating  technique  of 
kriging  “can  perform  well  in  image  enhancement  by  improving  the  resolution  of  satellite  imagery” 
and  “can  produce  images  which  are  more  faithful  than  cubic  convolution  to  the  original  image”  but 
leave  the  original  image  unchanged  (19:4-30). 

This  chapter  is  a  review  of  literature  dealing  with  the  basic  makeup  of  digital  imagery,  res¬ 
olution  of  imaging  systems,  current  digital  enhancement  techniques  and  the  application  of  kriging 
as  applied  to  digital  image  enhancement. 


2.2  Digital  Imagery 

Mather  (17:153)  defines  a  digital  image  as  “a  numerical  record  of  the  radiance  leaving  each  of 
a  large  number  of  rectangular  areas”  called  picture  elements  (pixels).  The  radiance  values  can  be 
stored  as  six,  seven,  eight,  or  ten-bit  words  for  later  access.  (17:153).  As  examples,  the  Advanced 
Very  High  Resolution  Radiometer  (AVHRR)  on  the  Television  and  Infrared  Observation  Satellite 
(TIROS)  weather  satellites  use  a  ten-bit  word  where  the  Landsat  Thematic  Mapper  uses  eight- 
bit  words.  All  pixels  discussed  for  particular  techniques  in  this  research  will  be  eight-bit  words. 
Associated  with  each  pixel  is  a  value  that  corresponds  to  the  relative  radiance  of  the  area  being 
studied.  The  pixel  typically  assumes  a  gray  (or  brightness)  level  value  ranging  from  0  to  255  where 
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0  represents  black  and  255  white.  A  digital  image  can  be  considered  a  three-dimensional  array 
where  the  x  and  y  axes  are  the  spatial  representations  of  the  image  and  the  z  axis  is  the  radiance 
value  of  the  pixel  (17:153).  A  computer  displays  the  monochromatic  image  on  an  output  device  by 
assigning  the  gray  level  to  pixels  on  the  device. 

2.2.1  Resolution  Enhancement.  Magnification  of  an  image  is  enlargement  of  part  of  an 
image  to  increase  its  resolution.  Resolution  is  the  degree  of  discernible  detail  and  is  dependent  on 
its  dimensions  (13:33).  The  goal  of  the  m.agnification  is  to  make  features  noticeable  that  were  not 
noticeable  in  the  original  image  (19:2-2).  A  simple  way  to  enhance  the  resolution  of  an  image  is 
to  copy  the  original  gray  level  pixel  values  and  duplicate  each  pixel  by  a  factor  called  the  zoom  or 
magnification  factor.  The  magnification  factors  are  usually  limited  to  two,  four,  and  eight  (17:95). 
This  makes  for  crude  images  that  can  be  worse  than  the  original.  A  better  technique  than  this 
copy-and-repeat  method  would  be  an  interpolator. 

2.3  Imaging  System  Resolution 

One  method  to  increase  the  spatial  resolution  of  an  image  is  to  improve  the  optical  system 
used  to  collect  the  EM  radiation  that  forms  the  image.  Larger  optics  and  advanced  CCD  arrays  will 
decrease  the  field  of  view  (FOV)  of  each  sensor  element  but  not  without  problems  and  limitations. 

2.3.1  Airy  Pattern.  The  irradiance  on  the  focal  plane  of  an  imaging  system  from  a  point 
source  takes  the  shape  of  a  Bessel  function  or  “Airy  pattern”  (named  after  Sir  George  Airy)  not 
a  point  as  one  might  expect  (10:137).  See  Figure  2.1  for  a  graphical  representation  of  an  Airy 
pattern.  This  diffraction  pattern  forms  a  three-dimensional  Airy  disk  (Figure  2.2),  on  the  sensor 
array  with  84  percent  of  the  EM  radiat'on  entering  the  sensor  contained  in  the  center  ring. 

The  radius  r  of  the  center  circle  depends  on  the  diameter  D  and  focal  length  /  of  the  collecting 
optics  and  the  wavelength  A  of  the  collected  radiation  by  the  following  relationships  (10:138): 

1.22/A/D  (2.1) 

which  can  be  represented  as: 

a  «  1.22A/D  (2.2) 

with  a  representing  the  angular  resolution. 
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Figure  2.1.  Airy  Pattern.  (10.139) 


2.3.S  Rayleigh  Criterion.  Spatial  resolution  can  then  be  defined  as  the  minimum  separation 
between  the  center  peaks  of  two  adjacent  Airy  disks.  This  “Rayleigh  criterion”  (Figure  2.3)  can  be 
applied  to  satellite  images  and  is  considered  the  best  possible  resolution  that  an  imaging  system 
can  attain.  A  satellite  at  an  altitude  A,  viewing  the  surface  of  the  earth  with  an  angle  4>  from 
nadir,  has  a  linear  separation  of  S  on  the  ground  given  by  (10:140): 

S«  (1.22/AA/D)sec(^  (2.3) 

2.3.2. 1  Diffraction' Limited.  The  Rayleigh  criterion  represents  the  best  possible  res¬ 
olution  an  imaging  system  can  achieve  and  is  called  diffraction-limited  resolution.  In  reality,  res¬ 
olution  will  be  limited  by  many  other  factors  besides  diffraction.  Physical  size  and  spacing  of 
the  detector  elements  in  an  array  place  additional  limits  on  the  imaging  system.  Elements  must 
be  small  enough  to  distinguish  between  the  two  Airy  disks.  Further  reductions  in  resolution  are 
caused  by  lack  of  stability  in  the  space  platform,  observation  geometry  providing  non-uniform  FOV, 
atmospheric  turbulence,  and  aberrations  in  the  focusing  mirrors. 

To  obtain  improvement  in  the  resolution  of  the  imaging  system,  alterations  in  focal  length, 
diameter  of  optics,  CCD  element  size,  satellite  altitude,  or  observed  wavelength  are  required.  Figure 
2.4  illustrates  the  geometry  of  the  FOV  of  one  pixel  in  a  satellite  imaging  system.  Since  minimum 
altitude  is  limited  by  orbital  mechanics  and  wavelength  is  controlled  by  the  specific  application  of 
the  system,  changes  in  the  optics  and  array  are  the  only  options  available  to  satellite  designers  to 
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Figure  2.2.  Airy  Disk.  (10:140) 


maximize  resolution.  Focal  length,  mirror  diameter,  or  CCD  element  size  may  be  altered  to  achieve 
the  required  resolution. 

2.3.3  Resolution  Example.  The  following  is  an  example  illustrating  the  limits  of  current 
state-of-the-art  imaging  systems.  This  example  is  adapted  from  class  notes  from  for  PHYS  621, 
Electro-Optical  Space  Systems  Technology  (11). 

2.3.3. 1  Element  Size.  If  the  Hubble  Telescope  were  pointed  toward  Earth  to  image 
and  distinguish  characters  on  an  automobile  license  plate,  what  size  detectors  and  optics  would  be 
required?  To  satisfy  the  Rayleigh  criterion,  the  system  must  resolve,  in  this  case,  objects  (plate 
characters)  2  cm  apart,  so  using  the  following  Hubble  optics  dimensions: 

D  =  2.4m 
/  =  6.4m 
A  =  700km 
S  =  2cm 

and  the  previous  equations  would  result  in: 

a  full  =  S/A 

=  .02m/7  X  10*m 


Figure  2.3.  Rayleigh  Criterion.  (10:140) 

=  2.86  X  10-®rad 


then 


r  =  q/ 

=  (2.86  X  10“*rad)(6.4tn) 

=  0.183/im 

These  results  show  that  the  CCD  elements  must  be  smaller  in  size  than  the  wavelength  of  the 
received  light  (visible  range  of  .4/im  to  .7/im)  in  order  to  distinguish  characters  on  the  plate. 
Current  state-of-the-art  Si  CCD  elements  are  13/im  in  size. 

2. 3. 3. 2  Diffraction-Limited.  For  diffraction- limited  resolutions  (the  best  resolution 
possible  for  the  Hubble  Telescope),  the  smallest  distinguishable  distance  projected  upon  the  CCD 
array  is: 


r  «  1.22/A/D 

«  1.22(6.4m)(.5^m)/2.4m 

«  1 .63#im 

(2.4) 
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Figure  2.4.  Imaging  System  Geometry.  (10:140) 


This  result  is  larger  than  the  r  necessary  to  identify  the  plate  numbers.  In  other  words,  the  Hubble 
Telescope’s  mirror  must  be  altered  to  perform  this  task.  To  accomplish  this  task,  the  Hubble 
Telescope’s  mirror  must  have  the  following  diameter  D: 

D  «  1.22/A/r 

«  1.22(6.4m)(.183^m)/2.4m 
sa  21.3m 

The  diameter  of  optics  to  perform  the  task  would  be  unmanageable  in  space  as  well  as  being  beyond 
current  technology. 

This  example  demonstrates  the  changes  required  to  alter  a  satellites  optical  system  to  increase 
the  resolution  of  the  imagery  produced  by  the  system  The  problems  associated  with  such  a  change 
include: 

1.  CCD  element  size.  State-of-the-art  CCD  technology  cannot  produce  the  necessary  element 
size. 

2.  Mirror  diameter.  To  compensate  for  the  larger  element  size,  the  mirror  of  the  satellite  must 
be  very  large  which  would  make  it  unmanageable  in  space  as  well  as  being  beyond  current 
technology. 
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2-4  Current  Resolution  Enhancement  Techniques 


An  alternative  to  altering  the  optics  of  an  imaging  system  is  to  use  ground  software  to  enhance 
the  resolution  of  an  image.  This  method  is  considerably  cheaper  and  is  without  the  problems 
associated  with  increasing  the  power  of  the  optics.  For  software  magnification,  the  portion  of 
interest  within  the  image  is  selected  and  this  area  is  expanded  (pixels  are  spaced  out)  and  new 
pixels  are  added  between  the  original  pixels.  Gray  levels  are  then  estimated  and  assigned  to  the 
coordinates  of  the  new  pixels.  Currently  three  processing  methods  are  used  today  to  estimate  pixel 
values  and  they  are: 


•  Nearest-neighbor 

•  Bilinear  interpolation 

•  Cubic  convolution 


2.4.1  Nearest- Neighbor.  The  nearest-neighbor  technique  of  image  magnification  is  the  sim¬ 
plest  method  in  use  today.  The  value  of  the  new  pixel  takes  the  value  of  the  original  pixel  nearest 
to  its  location.  It  assigns  the  gray  level  of  the  closest  pixel  (DNkj)  to  the  new  pixel  (Zu)  (3:734). 
The  nearest  neighbor  is  found  by: 


and  then 


k  =  integer  part  of  (i-f-O.S) 
I  —  integer  part  of  {j-l-0.5) 


Zij  =  Zki 


(2.5) 

(2.6) 


(2.7) 


Nearest-neighbor  hats  two  main  advantages:  it  is  fast  (takes  little  processor  time)  and  ensures 
that  the  pixel  values  in  the  enlarged  image  directly  relate  to  the  original  image  (17:138).  This 
technique  is  not  without  its  problems.  Bernstein  states  that  “the  major  drawback  to  nearest 
neighbor  assignment  is  the  discontinuities  which  are  introduced  by  its  zero-order  interpolation” 
which  causes  a  stairstep  or  blocky  appearance  to  the  image  (1:6-4). 


2.4.2  Bilinear  Interpolation.  The  bilinear  interpolation  method  uses  the  gray  levels  of  the 
four  nearest  neighbors  to  estimate  the  new  pixel  values.  Many  algorithms  to  implement  bilinear 
interpolation  are  possible  (1:6-4). 
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Bernstein  describes  an  algorithm  developed  by  IBM  which  minimizes  the  programming  com¬ 
plexity  and  maintains  high  accuracy; 

Zuv  —  Zii  -I-  d{Zi2  —  Zii)  +  d'[Z2i  +  d{Z22  —  Z21)  —  d{Zi2  —  Zn)]  (2-8) 

where  is  the  intensity  of  the  new  pixel  and  Zn,  Z\2,  Z21,  and  Z22  are  the  four  nearest 
neighbors  and  are  selected  as  shown  in  Figure  2.5.  This  technique  fills  the  desired  pixels  and 


9- Input  Inuge  umples 

X  -  Output  biuge  Sample  Mapped  Into  Image  Space 

Figure  2.5.  Bilinear  Interpolation.  (1:6-5) 

produces  smoother  images  but  tends  to  blur  the  image  (20:112). 

S.4.3  Cubic  Convolution.  Cubic  convolution  got  it  name  because  it  is  “based  on  the  fitting 
of  a  two-dimensional,  third-degree  polynomial  surface  to  the  region  surrounding  the  point”  (17:141). 
Gonzalez  says  “smoother  results  can  be  obtained  by  using  more  sophisticated  techniques,  . . .  which 
fit  a  surface  of  the  (sina;)/a;  type  through  a  much  larger  number  of  neighbors  . . .”  (13:300).  This 
means  that  the  16  nearest  pixel  gray  levels  are  used  to  estimate  the  value  of  the  new  pixel.  This 
technique  gives  a  more  natural-looking  image  without  the  blockiness  of  the  nearest-neighbor  or 
the  over-smoothing  of  the  bilinear  method  (17:141).  One  problem  in  the  implementation  of  the 
adaptation  of  this  method  is  that  it  requires  an  infinite  amount  of  data  to  accurately  estimate  the 
(sin  i)/*  function  (2:55). 
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Bernstein  describes  an  eilgorithm  used  for  cubic  convolution  which  uses  four  input  points 
(horizontal  interpolation  points)  and  16  neighbors  in  the  following  way; 

Z'k  —  ^*i[4  —  8(1  +  d)  +  5(1  +  d)^  —  (1  +  d)^] 

+■^*2(1  —  2d^)  +  d®) 

+Ml-2(l-d)2  +  (l-d)3] 

+Zt4[4  -  8(2  -  d)  +  5(2  -  df  -  (2  -  d)3] 

=  d{d[d(Zfc4  —  Zn;3  + Z|:2  —  Zfci) 

+{—Zk4  +  Ik3  —  2Zk2  +  2Zki)\ 

+(^lfc3  —  Zkl)}  +  Zk2 

where 

Z'k  =  d{d[d(Zk4  —  Zk3  +  Zk2  ~  Zki)  +  (Zk3  —  Zk4  (2-11) 

—2Zk2  +  2^j;i)]  +  {Zk3  —  Zki)}  —  Zk2 

then  the  final  pixel  values  are  solved  using  the  following  equation: 

Zuv  =  d'{d'[d'(^4  -  4  +  ^2  - (2.12) 

MZ'^  -  Z'  -  2Z'  +  2Z{)]  +  (Z'  -Z[))-  Z'^ 

See  Figure  2.6  for  an  illustration  of  the  pixel  configuration. 

2.5  Kriging 

Kriging  is  a  spatial  estimation  technique  developed  by  Matheron  in  the  early  1960s  to  estimate 
ore  reserves  from  spatial  core  samples.  Journel  describes  kriging  as  “a  local  estimation  technique 
which  provides  the  best  linear  unbiased  estimator  of  the  unknown  characteristics  studied”  (15:304). 
Clark  describes  kriging  as  “the  optimal  prediction  in  space  using  observations  taken  at  known  nearby 
locations”  and  is  a  weighted  average  where  the  weights  depend  on  the  location  of  the  points  used 
in  the  prediction  (7:240).  The  estimator  is  of  the  form: 

Z*  =  UiZi  +  W2^2  +  +  •  •  •  +  UfiZn  (2.13) 


(2.9) 


(2.10) 
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Figure  2.6.  Cubic  Convolution.  (1:6-7) 


where  the  w„  are  the  weights,  are  the  values  of  the  surrounding  data,  and  Z*  is  the  unbiased 
estimator  (5:99).  If  the  weights  used  for  the  estimates  sum  to  one  and  there  is  no  drift,  the  estimates 
are  considered  unbiased  (9:384).  The  weights  are  found  by  using  the  following  kriging  equations  in 
general  matrix  notation: 


0 

1 

1 

1 

X 

1 

1 

7(Aii) 

7(^12) 

■ ■ ■  7(*ln) 

Wi 

7(/>ip) 

1 

7(^21) 

7(/»22) 

•  •  •  7(fi2n) 

W2 

= 

7(/«2p) 

1 

7(/>„i) 

7(h„2) 

7{hnn) 

Wn 

s 

1 . 

(2.14) 
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where  7(hij)  are  the  calculated  semi-variograms  between  the  points  in  the  neighborhood,  y{hip) 
are  calculated  semi-variograms  between  the  points  in  the  neighborhood  and  the  estimated  point, 
and  A  is  a  Lagrange  multiplier  used  to  solve  the  simultaneous  equations. 

In  this  research,  kriging  is  used  for  spatial  prediction  of  image  pixel  values  during  the  en¬ 
hancement  of  imagery  resolution. 

2.5.1  Structural  Analysts.  Structural  analysis  attempts  to  find  the  size,  shape,  orientation, 
and  spatial  arrangement  of  samples  and  constitutes  the  variation  of  the  pixels  in  an  image  (9:239). 
Determination  of  the  variogram,  anistropy,  and  trend  constitute  the  structural  analysis  of  an  image. 

2. 5. 1.1  Variogram.  Any  discussion  of  kriging  must  start  with  a  discussion  of  the  var¬ 
iogram  which  represents  the  degree  of  continuity  or  correlation  of  a  data  set  (18:1250).  The  vari¬ 
ogram  describes  the  variance  of  the  difference  of  samples  within  the  data  set  and  is  calculated  by 
the  following  equation: 

i  (2.15) 

”  «=i 

where  z{x)  is  the  value  of  the  data  at  point  x  and  2(x-l-/>)  is  the  value  at  a  point  a  lag  h  from  x.  The 
semi-variogram  is  one-half  the  value  of  the  variogram.  We  now  have  a  measure  of  the  correlation 
between  data  values  h  apart.  As  will  be  seen  later,  these  calculations  are  easily  displayed  in 
graphical  form  with  lags  plotted  on  the  horizontal  axis  and  the  calculation  of  the  semi-variogram 
on  the  vertical  axis. 

Variogram  Example.  The  following  example  is  used  to  illustrate  the  calculation  of  an  experi¬ 
mental  semi-variogram  (adapted  from  Clark  (5:11-16)). 

The  data  shown  in  Figure  2.7  represents  a  stratiform  iron  deposit  in  which  holes  are  bored 
perpendicular  to  the  dip  (direction  ore  deposit  travels  into  the  earth)  of  the  ore  (see  Figure  2.8).  In 
this  case,  h  depends  on  the  distance  between  sample  pairs  and  the  orientation  of  drilling  provides 
a  convenient  square  grid  with  holes  100  feet  apart  (A).  The  data  given  is  the  average  value  of  the 
iron  (percentage  by  weight)  over  the  bore  sample. 

Semi-variogram  calculations  can  take  place  and  will  be  calculated  for  multiples  of  100  feet 
in  the  east-west  (0  degree)  direction.  When  h  —  0,  7*(0)  =  0.  When  h  =  100  feet,  all  pairs  100 
feet  apart  in  the  0  degree  direction  must  be  found  (Figure  2.9)  and  the  semi-variogram  7*  (100)  is 
calculated  by  Equation  2.16. 
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Figure  2.7.  Data  on  a  grid  for  the  calculation  of  an  experimental  semi-variogram.  (5:11) 


7*(100)  =  [(40  -  42f  +  (42  -  40)=*  +  (40  -  39)2  +  (39  -  37f  ^2.16) 

+(37  -  36)2  ^  (43  _  42)2  (42  _  39)2  4.  (39  _  39)2 

+(39  -  41)2  4.  (41  _  40)2  4.  (40  _  33)2  4.  (37  _  37)2 
+(37  -  37)2  4.  (37  _  35)2  4.  (35  _  33)2  4.  (33  _  37)2 
+(37  -  37)2  4.  (37  _  33)2  4.  (33  _  34)2  4.  (35  _  33)2 
+(35  -  37)2  4.  (37  _  30)2  4.  (35  _  35)2  4.  (35  _  35)2 
+(36  -  35)2  4.  (35  _  30)2  4.  (30  _  35)2  4.  (35  _  34)2 
+(34  -  33)2  4.  (33  _  32)2  4.  (32  _  39)2  4.  (29  -  28)2 
+(38  -  37)2  4.  (37  _  35)2  4.  (29  _  39)2 
+(30  -  32)2]  ^  (2  X  36) 

7*  (100)  =  1.46 

This  gives  a  value  for  7*  at  h  =  100  for  the  graph.  Now,  find  all  points  200  feet  apart  (Figure  2.10) 
and  calculate  7*  (200)  as  in  Equation  2.17. 


7*(200)  =  [(44  -  40)2  +  (40  -  40)2  +  (42  -  39)2  +  (40  -  37)2  (247) 
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Figure  2.8.  Cross-Section  Through  the  Iron  Ore  Deposit.  (5:12) 


-|-(39  -  36)2  ^  ^42  _  43)2  ^43  _  39^2  4.  ^42  _  39)2 

+(39  -  41)2  +  (39  -  40)2  4.  (4j  _  33^2  4.  (37  _  37^2 

+(37  -  35)2  4.  (37  _  33)2  4.  (35  _  37)2  4.  (33  _  37)2 

+(37  -  33)2  4.  (37  _  34)2  4.  (33  _  35)2  4.  (35  _  33)2 

+(37  -  36)2  4.  (36  _  35)2  4.  (35  _  35)2  4.  (35  _  35)2 

+(36  -  34)2  4.  (35  _  33)2  4.  (34  _  32)2  4.  (33  _  29)2 

+(32  -  28)2  4.  (33  _  35)2  4.  (35  _  39)2  4.  (39  _  29)2 

+(29  -  32)2]  ^  (2  X  33) 

T*(200)  =  3.30 

This  gives  a  second  point  to  graph.  Calculations  can  continue  up  to  /i  =  800  feet  but  Clark 
states,  “In  practice,  we  rarely  go  past  about  half  the  total  sample  extent.”  North-south  (90  degree) 
calculations  can  be  made  and  both  plotted  on  the  same  graph  (Figure  2.11). 

If  there  is  a  large  difference  between  the  plots  of  the  two  directions,  this  means  that  more 
information  is  needed  to  find  the  eixis  of  the  anisotropy  (this  will  be  described  later). 

Sill  And  Range.  As  samples  become  farther  apart,  the  correlation  between  these  points 
decrease  and  the  calculated  semi-variogram  value  will  increase.  Eventually,  in  the  ideal  case,  the 
samples  will  become  independent  and  the  semi-variogram  value  will  remain  constant  (5:6).  The  lag 
at  which  the  samples  become  independent  is  called  the  range  (a)  and  the  calculated  semi-variogram 
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Figure  2.9.  Identifying  all  Pairs  at  100  ft  Apart  in  the  East- West  Direction.  (5:13) 


where  this  occurs  is  called  the  sill  (C).  The  sill  is  where  the  correlation  between  data  points  is 
negligible  and  defined  as  the  ordinary  variance  of  the  sample  data  (5:7): 


1=1 


where 


i=l 

Nugget  Effect.  Davis  describes  the  nugget  effect: 


(2.18) 


(2.19) 


“. .  .the  nugget  effect  arises  because  the  regionalized  variable  is  so  erratic  over  a  very 
short  distance  that  the  semi-variogram  goes  from  zero  to  the  level  of  the  nugget  effect 
in  a  distance  less  than  the  sampling  interval.” (9:246) 

The  nugget  effect  (C'o)  is  random  in  nature  and  is  represented  by  a  spherical  model  with  a  very 
small  range  (5:9). 

Models. The  experimental  data  above  must  be  fit  to  a  model  to  estimate  points  for  the  later 
kriging  equations.  A  model  that  matches  the  experimental  data  must  be  selected  Jind  fit  to  the 
semi-variogram  curve. 
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Figure  2.10.  Identifying  all  Pairs  at  200  ft  Apart  in  the  East- West  Direction.  (5:14) 


Linear  Model.  The  data  for  the  previous  ore  example  suggests  that  a  linear  model  would 
fit.  The  linear  model  is  the  simplest  with  only  the  slope  (m)  describing  the  shape.  The  following 
equation  describes  its  shape. 

I  mh  -4-  Co  for  /»<  a 
7(/»)  =  < 

C  for  ft  >  a 

For  lags  less  than  the  range,  the  linear  model  makes  a  satisfactory  estimate  of  the  semi-variogram 
(9:246-247).  The  shape  of  the  linear  model  is  displayed  in  Figure  2.12. 

Spherical.  An  ideal  model  to  describe  the  semi-variogram  begins  at  the  origin  and  rises  to  a 
limit  (the  sill)  and  then  levels  out  to  a  constant  value  (9:246).  The  spherical  model,  described  by 
the  following  equation,  matches  this  description: 


j(h)  = 


I 


C[(3/2)(/i/a)  -  (l/2)(/»3/a3)]  -|-  Co  for  h  <  o 
C  -I-  Co  for  /»  >  o 


The  range  is  the  corresponding  lag  (ft)  where  the  sill  occurs.  Figure  2.13  represents  the  shape  of 
the  spherical  model  showing  the  range  and  sill. 
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Figure  2.11.  Experimental  Semi-Variograms  in  0®  and  90®  Directions.  (5:15) 


Figure  2.12.  Linear  Model. 
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Figure  2.13.  Spherical  Model.  (19:2-13) 


Exponential.  The  exponential  model  is  similar  in  nature  to  the  spherical  model  though  it 
never  levels  out  to  a  constant  value  or  sill  and  is  described  by  the  following  equation: 

7(A)  =  [1  -  6"*/“]  -I-  Co  (2.20) 

Figure  2.14  shows  the  exponential  model. 

2. 5. 1.2  Trend.  Trends  in  the  data  occur  when  samples  in  a  particular  direction  demon¬ 
strate  regional  fluctuations.  Since  stationarity  (the  mean  and  covariance  of  the  regionalized  variable 
are  known  to  exist  and  are  constant  over  the  area  (8:92))  is  assumed  prior  to  kriging,  the  trend 
must  be  removed  from  the  sample  data  to  prevent  biased  results  (5:120).  The  trend  can  be  removed 
by  fitting  a  linear  or  polynomial  function  to  the  data  and  performing  the  structural  analysis  on  the 
remaining  residuals. 

2. 5. 1.3  Anisotropy.  Anisotropy  explains  differences  in  semi-variograms  calculated  in 
different  directions  and  results  from  non-isotropic  data.  Geometric  anisotropy  occurs  when  vari¬ 
ations  in  semi-variogram  values  in  differing  directions  are  the  same  for  lags  that  are  multiples  of 
each  other  (19:2-15).  The  semi-variograms  will  have  the  same  sill  but  differing  ranges.  A  correction 
factor,  k,  can  be  used  to  equalize  the  ranges.  The  spherical  model  can  then  be  corrected  as  in  the 
following: 
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Figure  2.14.  Exponential  Model. 


to 

=  c(^-^)  +  C,  (2.22) 

which  makes  the  two  equations  equal. 

Figure  2.15  graphically  demonstrates  geometric  anisotropy  by  plotting  zero  and  ninety-degree 
semi-variograms  for  an  image  from  Channel  One  of  the  SPOT  multi-spectral  sensor.  The  differences 
between  the  two  plots  is  due  to  the  anisotropy  and  represents  the  non-isotropic  nature  of  the  image. 

2.5.2  Ordinary  Kriging.  Ordinary  kriging,  a  form  of  punctual  kriging  (predicting  values  of 
a  stationary  regionalized  point  variable),  solves  for  the  new  spatial  data  by  weighting  the  nearest 
neighbors  (size  of  16,  20,  25  . . .  dependent  on  the  range  of  influence)  and  substituting  into  the 
estimator.  Ordinary  kriging  solves  for  the  pixel  values  by  the  following  estimator  equation: 

n 

Z*  =  (2.23) 

i=l 
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SPOT  Channel  One  Image 


Figure  2.15.  Experimental  Semi-Variograms  in  0®  and  90®. 


where 

n 

=  1  (2.24) 

<=i 

If  the  weights  sum  to  one  and  no  trend  exists,  then  Z*  is  considered  an  unbiased  estimator.  This  is 
a  linear  estimator  because  it  is  a  linear  combination  of  the  samples  (5:99).  To  determine  the  weights 
necessary,  a  system  of  equations  must  be  derived.  If  there  are  n  samples  within  the  neighborhood 
used  to  estimate  a  point,  the  following  linear  equations  will  result: 


0 

-1- 

Wi 

+ 

W2 

+ 

W3 

-1- 

...  -1- 

W„  =  1 

A 

+ 

wa{hu) 

+ 

t«27(Al2) 

+ 

«'37(Ai3) 

-1- 

...  + 

Wn7(^ln)  =  7(*lp) 

A 

+ 

wi'r(h2\) 

+ 

lU27(/»22) 

+ 

tV3y(h23) 

+ 

...  q. 

W'n7(/»2r.)  =  7(h2p) 

A 

+ 

waihai) 

+ 

«^27(/»32) 

+ 

«'37(f‘33) 

+ 

...  q. 

=  7(/»3p) 

•f 

+ 

+ 

+ 

+ 

:  =  : 

A 

+ 

wnihni) 

+ 

«'27(Af.2) 

U)37(  ftns) 

-f 

...  + 

U^n7(^nn)  —  7(^np) 

where  is  the  weight  at  point  i  of  the  neighborhood  to  the  estimated  point  p  and  7(/»ij)  is  the 
semi-variogram  estimate  from  points  i  to  j.  The  regular  grid  nature  of  the  images  simplify  the 
calculations  of  the  semi-variograms  limiting  the  number  of  equation  systems.  These  equations  are 
solved  for  the  weights  and  used  to  estimate  the  new  pixel  value  Z*  as  in  Equation  2.23. 
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To  solve  the  system  of  equations,  matrices  are  developed  as: 


0  1  1 

1  7(/>ii)  7(*12) 

1  7(/»2i)  7(A22) 

1  7(/»r.l)  7(/ln2) 


1 

A 

1 

7(/>ln) 

uq 

7(/»ip) 

7(A2n) 

W2 

= 

7(/«2p) 

7(^nn) 

Wn 

7(/*np) 

(2.25) 


The  weights  are  calculated  for  all  estimated  points  and  equations  are  generated  to  produce  the 
estimator. 


2.5.3  Ordinary  Kriging  Example.  Davis  demonstrates  the  technique  of  ordinary  kriging  by 
estimating  the  elevation  of  the  water  table  at  Point  p  on  the  map  in  F'igure  2.16  (9:386-389).  In 
estimating  Point  p,  three  measured  elevations  will  be  used,  with  their  coordinates  given  in  Table 
2.1.  The  semi-variogram  produced  from  this  data  is  linear  with  a  slope  of  4.0  m^/km  for  a 
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Figure  2.16.  Water  table  elevations  (meters)  at  three  observations  wells.  (9:394) 


neighborhood  of  20  km.  The  distances  (/i)  between  the  points  are  shown  in  Table  2.2  with  the 
corresponding  semi-variograms  values,  7(Ai;),  listed  in  Table  2.3.  The  kriging  matrix  will  be  as 
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Table  2.1.  Water  Table  Elevation  Data  for  the  Ordinary  Kriging  Example. 


Location 

X  Coordinate 

Y  Coordinate 

Water  Table 
Elevation 

■l^n 

3.0 

4.0 

6.3 

3.4 

2.0 

1.3 

Point  p 

3.0 

3.0 

Table  2.2.  Distances  Between  Wells  and  Point  p  (Ordinary  Kriging). 


Location 

Well 

1 

Well 

2 

Well 

3 

Point 

P 

HZ  y 

0 

3.35 

2.88 

y 

0 

4.79 

3.32 

Q 

0 

1.97 

Table  2.3.  Semi-Variograms  for  Distances  Between  Wells  and  Point  p  for  the  Ordinary  Kriging 
Example. 


Location 

Well 

1 

Well 

2 

Well 

3 

Point 

P 

Well  1 

0 

13.42 

11.52 

4.00 

Well  2 

0 

19.14 

13.30 

Well  3 

0 

7.89 
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follows: 


0.0  1  1  1 

A 

1 

1  7(/»ii)  7(^12)  7(^13) 

Wi 

7(/»ip) 

1  7(^21)  7(^22)  7(^23) 

W2 

7(/»2p) 

1  7(/>3i)  7(^32)  7(/*33) 

W3 

.  7(/»3p) 

Substituting  data  from  Table  2.3  into  the  previous  kriging  equations  results  in  the  following  system 
of  equations  in  matrix  form: 


0.0  1  1  1 

A 

1 

1  0.0  13.4  11.5 

IWl 

4.0 

1  13.4  0.0  19.1 

W2 

13.3 

1  11.5  19.1  0.0 

U>3 

7.9 

The  diagonal  of  the  matrix  has  all  zero  values  and  the  matrix  is  symmetrical  as  discussed  earlier, 
resulting  from  semi-variogram  calculations  7(h<j)  within  the  neighborhood.  Solving  this  system  of 
equations  result  in  the  following  weights: 


A 

-0.7267 

Wl 

0.6039 

W2 

0.0868 

W3 

0.3093 

To  predict  the  value  of  Point  p,  the  preceding  weights  are  substituted  into  the  estimator  and  the 
result  is  as  follows: 

3 

Zp  =  ^  ^  Wi  Zi 

i  =  l 

=  0.6039(120.0)  +  0.0868(103.0)  +  0.3093(142.0) 

=  125.3  meters 

Ordinary  kriging  can  be  used  if  no  drift  exists  within  the  image  or  sample  data. 

2.5.4  Universal  Kriging.  As  discussed  earlier,  the  data  must  be  stationary  (no  drift)  for  the 
kriging  process  to  be  effective  and  unbiased.  Davis  states: 
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“The  drift  consists  of  the  average  or  expected  vedue  of  the  regionalized  variable  within  a 
neighborhood,  and  is  a  slowly  varying,  nonstationary  part  of  the  surface.  The  residual 
is  the  difference  between  the  actual  measurements  and  the  drift.”  (9:393) 


After  the  drift  is  removed,  the  residuals  must  be  stationeiry  and  can  be  kriged.  After  kriging,  the 
residuals  are  combined  with  the  drift  to  retrieve  the  necessary  data.  Universal  kriging  is  similar 
to  ordinary  kriging  but  relates  the  semi-variogram  of  the  residuals  of  the  drift,  the  size  of  the 
neighborhood,  and  the  drift  model  into  the  equations.  From  this  interrelationship,  the  following 
matrix  results: 
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where  A,  and  V)  are  coordinates  of  the  neighborhood  points  with  Xp  and  Yp  the  coordinates 
of  the  estimated  point.  These  coordinates  are  measured  relative  to  a  control  point  within  the 
neighborhood. 


2.5.5  Universal  Kriging  Example.  Davis  demonstrates  the  technique  of  universal  kriging  by 
continuing  the  estimation  of  the  previous  example  for  the  elevation  of  the  water  table  at  point  p 
on  the  map  in  Figure  2.17  (9:386-389).  In  estimating  Point  p  for  this  case,  five  known  elevations 
will  be  used,  with  their  coordinates  given  in  Table  2.4. 


Table  2.4.  Water  Table  Elevation  Data  for  the  Universal  Kriging  Example. 


Location 

X  Coordinate 

Y  Coordinate 

Water  Table 
Elevation 

n 

3.0 

4.0 

6.3 

3.4 

2.0 

1.3 

Well  2A 

3.8 

2.4 

115.0 

Well  5 

1.0 

3.0 

148.0 

Point  p 

3.0 

3.0 
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Figure  2.17.  Water  table  elevations  (meters)  at  five  observation  wells.  (9:394) 


The  semi-variogram  produced  from  this  data  is  the  same  as  the  ordinary  kriging  and  is  linear 
with  a  slope  of  4.0  m^/km  for  a  neighborhood  of  20  km.  The  distances  (hjj)  between  the  points 
are  shown  in  Table  2.5  with  the  corresponding  semi-variogram  values,  y{hij),  listed  in  Table  2.6. 


Table  2.5.  Distances  Between  Wells  and  p  for  the  Universal  Kriging  Example. 


Location 

Well 

1 

Well 

2 

Well 

3 

Well 

2A 

Well 

5 

Point 

P 

0 

3.35 

2.88 

1.79 

2.24 

1.00 

Well  2 

0 

4.79 

2.69 

5.32 

3.32 

Well  3 

0 

2.11 

1.97 

1.97 

Well  2A 

0 

2.86 

1.00 

Well  5 

0 

2.00 
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Table  2.6. 


Semi-Variograms  for  Distances  Between  Wells  eind  Point  p  for  the  Universal  Kriging 
Example. 


Location 

Well 

1 

Well 

2 

Well 

3 

Well 

2A 

Well 

5 

Point 

P 

Well  1 

0 

13.42 

11.52 

7.16 

8.96 

4.00 

Well  2 

0 

19.16 

10.76 

21.28 

13.28 

Well  3 

0 

8.44 

7.88 

7.88 

Well  2A 

0 

11.44 

4.00 

Wells 

0 

8.00 

The  kriging  matrix  will  be  as  follows: 
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Substituting  data  from  Table  2.6  into  the  previous  kriging  equations  result  in  the  following  system 
of  equations  in  matrix  form: 
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Solving  this  system  of  equations  result  in  the  following  weights: 


A 

-0.7275 

0.0826 

02 

-0.2555 

Wi 

0.3797 

W2 

-0.0171 

t03 

0.1115 

W4 

0.4356 

Ws 

0.0903 

To  predict  the  value  of  Point  p,  the  preceding  weights  are  substituted  into  the  estimator  eind  the 
result  is  as  follows: 

Z*  =  0.3797(120)  -  0.0171(103)  +  0.1117(142)  +  0.4356(115)  +  0.0903(148) 

=  123.09  meters 


2.6  Conclusions 

This  chapter  reviewed  basic  concepts  in  resolution,  current  resolution  enhancement  tech¬ 
niques,  and  kriging.  Resolution  depends  upon  the  optics  of  an  imaging  system.  In  order  to  increase 
the  resolution  of  an  image,  the  optics  of  the  systems  must  be  adtered.  Extensive  alterations  to 
optical  and  sensing  systems  are  extremely  expensive  and  difficult  from  a  technological  standpoint. 

Ground  processing  of  images  can  improve  resolution  without  the  costs  associated  with  chang¬ 
ing  an  imaging  system.  Current  techniques  used  are  nearest-neighbor,  bilinear  interpolation,  and 
cubic  convolution.  Kriging  can  be  applied  to  enhance  resolution  of  digital  satellite  imagery  because 
it  is  an  unbiased  spatial  predictor.  Ordinary  and  universal  kriging  were  discussed  as  the  preferred 
methods  but  if  trends  exist  in  the  data,  universal  kriging  is  the  single  preferred  method. 
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III.  Methodology 


The  following  chapter  summarizes  the  steps  taken  to  complete  the  objectives  listed  in  previous 
chapters.  The  methodology  includes  data,  implementation  of  current  methods,  kriging  implemen¬ 
tation,  program  operation,  and  output. 

3.1  Image  Data 

The  image  data  used  in  this  research  was  collected  from  various  sources  ranging  from  satellites 
to  digitized  aerial  photographs.  The  predominant  image  type  used  was  from  the  French  Systems 
Pour  I’Observation  de  la  Terre  (SPOT)  satellite  system.  The  SPOT-1  satellite  has  a  circular,  sun- 
synchronous  orbit  with  an  altitude  of  832  km  and  an  inclination  of  98.7®  with  its  orbit  repeating 
every  26  days  (16:581-582). 

The  sensor  provides  multi-spectral  data  from  three  3000-element  arrays  at  20m  resolution 
digitized  in  eight-bit  words  and  transmitted  at  25  Mbps  (16:584).  Channel  One  views  over  a 
spectral  range  of  0.50  to  0.59/im,  Channel  Two  over  0.61  to  0.68/im,  and  Channel  Three  over  0.79 
to  0.89^m.  The  system  uses  a  push-broom  scanner  which  samples  the  detectors  along  the  array  and 
then  forms  the  image  by  assembling  successive  scans  as  the  satellite  travels  over  the  earth.  This 
system  eliminates  geometric  errors  due  to  variations  in  speed  caused  by  similar  mirror  systems  as 
well  improving  the  signal-to-noise  performance  by  increasing  the  dwell  over  the  target  area  (16:582). 
The  images  used  were  corrected  for  geometric  distortions  due  to  earth  curvature  and  nonlinearities 
in  the  sensor  sweep  (16:612). 

3.1.1  Image  Used.  Images  were  chosen  to  provide  a  variety  of  terrain  and  spectral  qualities. 
The  goal  was  to  test  the  techniques  with  as  many  varying  spectral  qualities  to  verify  kriging’s  and 
cubic  convolution’s  ability  to  handle  varying  imagery. 

Multi-spectral  data  was  used  from  Channels  One,  Two,  and  Three  from  SPOT’s  multi-spectral 
scanning  system.  Images  tested  were  taken  from  ail  channels  to  verify  each  techniques  ability 
to  handle  variations  in  spectral  information.  Sections  of  water  were  chosen  for  their  isotropic 
characteristics  and  urban  areas  were  selected  due  to  their  non-isotropic  nature.  This  was  done  to 
fully  test  kriging’s  ability  to  handle  many  image  characteristics.  A  cross  reference  of  image,  sizes, 
statistics  on  the  pixels  within  the  image  sections,  see  Appendix  5.2. 

3.1.2  Sub-Sampling.  The  image  sections  used  were  sub-sampled  to  test  the  techniques.  For 
the  purposes  of  this  research,  sub-sampling  consist  of  the  removal  of  every  other  row  and  column. 
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These  will  be  reconstructed  by  the  estimation  techniques  tested.  The  original  image,  four  times 
the  size  of  that  specified  when  running  the  program,  is  read  out  to  a  file  for  future  comparison. 


3. 2  Current  Methods 

Cubic  convolution  was  used  as  the  method  to  test  against  kriging. 

3.2.1  Neighborhood.  For  testing  purposes,  a  fixed  4  by  4  neighborhood  was  chosen  for  its 
simplicity  as  well  as  matching  that  neighborhood  used  by  the  cubic  convolution  algorithm.  In  an 
operational  program,  the  neighborhood  should  be  responsive  to  the  range  where  the  diagonal  radius 
will  equal  the  range  (19:3-10). 


3.3  Kriging  Implementation 

Both,  ordinary  and  universal  kriging  were  implemented  within  the  test  program  to  see  if  any 
differences  exist  between  the  two  techniques  for  image  resolution  enhancement.  A  flag  within  the 
program  can  be  set  to  change  kriging  method. 

3.3.1  Structural  Analysis.  Structural  analysis  describes  the  size,  shape,  orientation,  and 
spatial  arrangement  of  samples  and  constitutes  the  variation  of  the  pixels  in  an  image  (9:239). 
Usually,  this  relates  to  the  physical  description  of  an  ore  structure  underground  but  can  be  applied 
to  the  data  of  an  image.  Determination  of  the  variogram,  anisotropy,  and  trend  constitute  the 
structural  analysis  of  data.  Structural  analysis  was  performed  only  on  the  specific  image  sections, 
not  the  overall  image,  to  limit  the  trend  and  anisotropy  of  the  sub-sampled  image  and  maximize 
the  accuracy  of  the  final  results. 

3. 3. 1.1  Trend.  An  option  to  remove  trend  within  the  data  by  the  use  of  least  squares 
regression  method  was  included.  The  result  of  the  technique  was  a  set  residuals,  one  for  each 
data  point.  It  is  this  residual  data  that  creates  the  semi- variogram  and  subsequently  kriged.  The 
resulting  kriged  image  is  then  recreated  by  calculating  the  final  pixel  values  using  the  regression 
equation  and  the  residual  data. 

3. 3. 1.2  Anisotropy.  Zonal  and  geometric  anisotropy  was  no  corrected  for  in  this  re¬ 
search.  Preliminary  results  indicated  that  little  geometric  anisotropy  exist  within  the  images.  See 
Chapter  4  for  further  details  of  this  data. 
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3.3. 1.3  Nugget.  From  previous  research,  the  nugget  was  set  to  zero  for  all  calculations. 
McGee  states,  “a  model  which  fit  well  near  the  origin  of  the  variogram  [was  needed]  . . .  the  nugget 
effect  of  the  spherical  model  was,  therefore,  set  to  zero  (19:4-9).” 

3.3.2  Kriging  Equations.  Both  ordinary  and  universal  kriging  were  used  for  this  research 
to  determine  the  best  suited  method.  Both  techniques  are  used  to  solve  for  the  following  weight 
estimator: 


Z*  —  uiZi  -b  U2Z2  +  uaZs  -!-•••  +  UInZn 


which  will  be  used  throughout  the  image  to  calculate  the  necessary  pixel  values.  If  the  weights  sum 
to  one  and  no  trend  exists,  then  Z*  is  considered  on  unbiased  estimator. 

3.3.2. 1  Ordinary  Kriging.  If  no  trend  exists  within  the  image  or  section  of  image  to 
be  kriged,  the  simpler  form  of  kriging  called  ordinary  kriging  may  be  used  which  solves  for  the 
weights  using  the  system  of  equations  found  in  Section  2.5.2.  The  weights  are  calculated  for  all 
points  within  the  neighborhood  and  equations  are  generated  to  produce  the  estimator. 

3. 3. 2.2  Universal  Kriging.  If  data  has  a  drift,  universal  kriging  is  used  to  compensate. 
Universal  kriging  is  similar  to  ordinary  kriging  but  relates  the  semi- variogram  of  the  residuals  of 
the  drift,  the  size  of  the  neighborhood,  and  the  drift  model  into  the  equations.  In  this  research, 
pixel  one  in  the  upper  left  hand  side  (Figure  3.1)  represents  the  zero  control  point. 

3.4  Program  Operation 

See  Appendix  A  for  a  detailed  treatment  of  the  makeup  eind  operation  of  the  test  program. 

3. 4-0. 3  Weights.  Due  to  the  regular  grid  nature  of  the  images,  only  three  systems  of 
weights  were  needed  to  estimate  the  points.  These  weights  can  then  be  used  throughout  the  image 
section  to  and  estimate  all  necessary  points.  See  Figure  3.1  for  the  position  of  the  points  estimated 
within  the  neighborhood.  Due  to  the  size  of  the  neighborhood,  a  bordc'r  around  the  test  image  is 
required  to  estimate  all  points.  This  research  developed  a  border  within  the  test  image  to  supply 
these  points. 
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Figure  3.1.  Neighborhood  Pixels. 


3.5  Tests 

Comparison  tests  were  made  to  compare  universal  and  ordinary  kriging  with  regression,  krig- 
ing  with  cubic  convolution,  interpolation  with  the  spherical  and  exponential  model,  nugget  and  no 
nugget,  and  regression  and  no  regression. 

3.6  Statistical  Analysis 

To  compare  images,  statistical  test  are  required.  The  estimated  pixels  are  subtracted  from 
the  original  image  pixels  and  a  difference  is  found.  The  absolute  value  of  this  difference  is  found 
and  a  sample  mean  and  standard  deviation  of  the  differences  is  determined.  These  are  found  by 
the  following: 


El"*  I 

»=i 

n 
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and 


s  = 


^1 


*  =  1 


1  =  1 


n(n-  1) 


S.6.]  Variance  Test.  Variances  are  then  compared  to  see  if  they  are  considered  equal  to 
determine  which  t  test  to  use.  An  F  statistic  is  calculated  as  a  comparison  with  the  theoretical. 
The  experimental  statistic  is  calculated  by  the  following: 


If  Fo  >  you  reject  the  hypothesis  that  the  variances  are  equal.  This  is  the  critical 

region  for  testing  the  null  hypothesis  s{  =  sj  against  the  alternative  hypothesis  Sj  >  Sj.  To 
Lest  the  null  hypothesis  against  the  alternative  hypothesis  that  Sj  <  s?.  Ihen  Fg  <  where 

(21:263). 


3.6.S  Mean  Test.  With  the  means  determined,  two-sided  tests  concerning  the  mean  are 
constructed  (21:11).  A  t  statistic  is  calculated  to  perform  the  test  as  in  the  following  (12:241): 


t  = 


(^1  -  h) 


(3.1) 


We  assumed  that  we  are  sampling  normal  populations  and  that  we  can  apply  the  central  limit 
theorem  so  we  can  approximate  ci  and  <T2  with  si  and  s^  when  ni  and  ^2  are  both  greater  or  equal 
to  30  (12:242).  In  this  case,  they  are  equal  to  867  for  our  40  by  40  images. 


Hg  -.  ftl  =  H2 

(3.2) 

(3.3) 

The  results  of  these  tests  are  found  in  Appendix  C.5. 

Comparisons  are  made  between  cubic  convolution  and  kriging,  universal  kriging  and  ordinary 
kriging,  the  exponential  and  spherical  models  and  interpolation,  and  nugget/non-nugget  models. 
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3. 7  Summary 

Image  data  used  for  enhancement  was  sub-sampled  from  various  sources  including  SPOT 
imagery  and  digitized  aerial  photographs.  The  test  program  allows  the  operator  to  set  various 
flags  within  the  code  to  vary  the  test  conditions.  Some  of  the  conditions  that  can  be  changed 
are  the  semi-variogram  model  to  be  used  and  the  method  of  kriging  to  test.  Drift  is  removed 
by  least-squares  regression  technique  and  the  resulting  residuals  generate  semi-variograms  and  are 
then  kriged.  The  kriged  image  is  recombined  with  the  regression  equation  and  compared  with  the 
original  image.  Statistics  are  generated  on  the  comparison  and  are  compared  with  like  statistics 
from  the  cubic  convolution  technique  to  determine  the  best  technique. 
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IV.  Results 


Cubic  convolution,  ordinary,  and  universal  kriging  were  applied  to  various  digital  image  sec¬ 
tions  from  two  SPOT  satellite  images  each  with  20m  resolution  from  the  system’s  multi-spectral 
scanner.  See  Appendix  5.2  for  description  of  the  images.  Each  image  section  was  cut  from  two 
different  channels  with  the  following  spectral  ranges: 

•  Channel  2  -  0.61  to  0.68/im 

•  Channel  3  -  0.79  to  0.89/im 

Each  technique  was  used  on  all  channels  to  test  their  ability  to  handle  varying  spectral  conditions 
(since  each  spectral  channel  varies  significantly).  This  chapter  contains  results  from  initial  structural 
analysis  to  determine  the  proper  model  to  use  and  see  if  any  anisotropy  exists  within  the  images. 

4-1  Siructural  Analysis. 

Semi-variograms  of  residual  data  were  plotted  to  perform  a  structural  analysis  on  the  images 
to  determine  the  proper  model  to  incorporate  within  the  program  and  to  determine  whether  any 
anisotropy  exists  within  the  image.  Semi-variograms  are  displayed  in  Appendix  B.8. 

4.1.1  Model  Determination.  Most  semi-variograms  plotted  for  this  research  indicate  the  use 
of  no  specific  model  for  use  in  estimation  of  the  kriging  equations.  The  plots  show  linear  tendencies 
for  lags  below  the  m£Lximum  distance  of  a  4  x  4  neighborhood  (\/l8). 

4.1.2  Anisotropy.  The  differences  between  the  0®  and  90®  directions  are  minor,  especially 
for  lags  less  than  that  used  for  neighborhood  calculations,  suggesting  no  significant  anisotropy  exists 
for  the  images  tested. 

4.1.3  Range.  Most  ranges  for  the  semi-variograms  were  ideal  for  the  4x4  neighborhoods 
used  with  some  smaller  indicating  a  small  neighborhood  could  be  used  in  some  instances. 

4.2  Enhanced  Images 

With  the  modified  mean  calculation,  runs  were  repeated  for  the  image  previously  enhanced 
as  well  as  runs  on  the  remaining  image  sections  not  enhanced  in  the  initial  runs.  Image  pixel 
mean  and  standard  deviation  was  collected  from  the  original  image  data.  This  data  is  displayed  in 
Appendix  5.2. 
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4-2.1  Test  Statistics.  Tests  were  performed  to  test  universal  and  ordinary  kriging  as  well  as 
kriging  with  cubic  convolution. 

4.2.2  Comparisons  of  Two  Population  Variances.  To  determine  the  test  to  perform  on  the 
means  of  the  two  populations,  a  test  of  the  variances  was  performed  for  those  cases  to  be  compared. 
The  Fo  statistics  are  recorded  in  Appendix  C.5.  We  would  reject  Hg  (s?  =  s\)  against  (sj  >  Sj) 
if  Fg  >  f^a/2,t»i-i,n3-i-  ^a/2,ni-i,na-i  for  an  o  =  .05,  equals  1.16.  This  number  is  the  theoretical 
value  and  can  be  obtained  from  tables  or  a  statistical  computer  package.  We  reject  the  null 
hypothesis  that  the  variances  are  considered  equal  if  the  experimental  statistic  Fg  >  1.16.  For  the 
comparisons  of  krigi.ag  and  cubic  convolution,  for  all  but  two  variances,  the  hypothesis  is  rejected. 
For  the  comparison  of  interpolation  and  the  spherical  model,  10  reject  the  hypothesis  while  12 
accept  the  hypothesis.  For  the  comparison  of  interpolation  and  the  exponential  model,  11  reject 
the  hypothesis  while  11  accept  the  hypothesis. 

4-2.3  Comparisons  of  Two  Population  Means.  At  this  point  a  test  was  needed  to  compare 
the  means  to  see  if  they  were  statistically  different.  An  experimental  t  value  was  calculated  using 
the  experimental  means  and  variances.  A  95  percent  confidence  level  was  used  making  a  =  .05. 
The  two  sided  test  used  test  the  following  hypothesis: 


=  M2 

(4.1) 


against 


Ml  M2 


with  a  95  percent  confidence  interval  of  -<a/2,n-i  <t  <  <0/2, n-i-  If  t  is  within  the  interval,  you 
accept  Hg  and  if  not,  you  reject  Hg  and  accept  Hg.  Data  from  all  tests  are  displayed  in  Appendix 
C.5.  The  confidence  interval  for  the  data  is  -- 1.965  <  t  <  1.965.  If  the  test  statistic  (t)  is  within 
the  interval,  the  means  are  considered  equal  and  outside  the  interval,  they  are  not  equal. 

4.2.3. 1  Universal/Ordinary  Kriging  Comparison.  Runs  were  made  to  collect  data  on 
both  universal  and  ordinary  kriging.  Data  on  mean  difference  and  standard  deviation  were  collected 
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at  the  end  of  each  run.  No  difference  between  the  two  methods  was  noticed  to  two  decimal  points 
of  accuracy.  This  data  is  not  included  within  this  document. 

4. 2. 3. 2  Kriging/ Cubic  Comparison.  The  variances  for  kriging  were  consistently  smaller 
than  that  produced  by  cubic  convolution.  Many  were  not  statisticcdly  different  than  that  produced 
by  cubic  convolution  but  many  were  when  making  the  Fg  statistic  comparison. 

Data  on  mean  difference  and  standard  deviation  were  collected  at  the  end  of  each  run.  This 
data  is  displayed  in  Appendix  C.5.  Kriging  consistently  produced  means  better  than  cubic  convolu¬ 
tion  for  the  images  tested.  All  images  produced  by  kriging  had  lower  means  than  cubic  convolution 
though  only  four  of  the  images  produced  t  statistics  that  concluded  that  kriging  is  better.  All 
remaining  t  statistics  were  within  the  confidence  interval  concluding  that  the  means  were  the  same. 

4. 2. 3. 3  Inierpolation/Spherical  Model  Comparison.  Data  on  mean  difference  and  stan¬ 
dard  deviation  were  collected  at  the  end  of  each  run.  This  data  is  displayed  in  Appendix  C.5.  The 
data  indicates  that  there  is  no  significant  difference  between  the  techniques  with  the  interpolation 
method  producing  better  results  for  images  with  little  or  no  anisotropy. 

4. 2. 3. 4  Interpolation/ Exponential  Model  Comparison.  Data  on  mean  difference  and 
standard  deviation  were  collected  at  the  end  of  each  run.  This  data  is  displayed  in  Appendix  C.5. 
Results  of  these  test  were  similar  to  the  comparison  of  interpolation  and  the  spherical  model. 

4. 2. 3. 5  Nugget/No  Nugget  Comparison.  Images  produced  using  the  nugget  effect  have 
significantly  high  means  differences  than  that  produced  with  no  nugget  effect. 

4.2.4  Regression/No  Regression  Comparison.  Images  were  produced  to  see  if  'he  inclusion 
of  the  regression  section  of  the  code  was  necessary  to  improve  the  mean  difference.  Data  showed 
no  difference  between  the  kriging  with  regression  and  that  without.  This  data  is  not  included  in 
this  document  to  avoid  redundancy. 

4.2.4. 1  Regression/Universal  Kriging  Comparison.  Images  produced  by  the  elimina¬ 
tion  of  the  regression  function  and  then  using  universal  kriging  produced  similar  results  than  that 
produced  using  the  regression.  Only  slight  differences  (one  to  two  hundredths)  were  noticed.  This 
data  was  not  included  since  it  varied  little  with  the  original  data. 
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4-3  Summary 

Semi-variograms  produced  from  the  kriged  images  show  large  fluctuations  in  the  range  of 
influence  causing  difSculty  automating  the  structural  euialysis  function  of  the  program  though  little 
anisotropy  exists  within  the  images  especially  for  lower  lags.  The  lower  lags  are  used  to  estimate 
the  kriging  equations.  No  differences  between  the  use  of  universal  and  ordinary  kriging  were  seen 
when  drift  was  removed  prior  to  kriging.  Kriging  consistently  produced  images  with  lower  means 
differences  than  that  produced  by  cubic  convolution.  Twenty-seven  percent  of  the  images  tested 
produced  significantly  better  meeins  for  kriging  while  45  percent  of  the  images  produced  lower 
variances.  Comparisons  did  not  produce  significant  differences  between  the  different  models  while 
the  inclusion  of  a  nugget  produced  larger  mean  differences.  The  final  comparison  between  the  use  of 
regression  and  without  it  (using  universal  kriging)  produced  images  with  the  same  means  as  those 
produced  with  the  regression  function  included  in  the  program.  This  would  speed  the  processing 
considerably  since  the  regression  is  the  most  time  consuming  aspect  of  the  program. 
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V.  Recommendations 


The  recommendations  in  this  chapter  suggest  improvements  to  the  image  processing  program 
to  further  improve  the  speed  and  performance  of  the  dedicated  kriging  image  enhancement  as  well 
as  a  suggestion  for  a  further  use  of  the  semi-variogram  plots. 


5.1  Program  Improvements. 

Several  improvements  can  be  made  to  the  program  that  may  improve  kriging’s  ability  to 
enhance  resolution  of  digital  images  as  well  as  improve  the  speed  of  the  processing. 

5.1.1  Model  Estimation.  The  spherical  model,  exponential  model,  and  interpolation  meth¬ 
ods  were  used  to  estimate  the  semi-variogram  values  within  the  kriging  equations.  The  semi- 
variogram  plots  show  for  the  lags  of  that  used  to  make  the  estimations,  a  linear  model  could  be 
used.  This  is  why  the  linear  interpolation  was  attempted.  A  linear  model  could  be  incorporated 
into  the  program  to  aid  in  this  process. 

5.2  Summary 

Several  things  could  be  done  to  improve  the  speed  of  the  processing  as  well  as  improving  the 
results  obtained  from  kriging.  Faster  processors  and  different  software  can  be  used  to  increase  the 
speed  and  handle  the  data  better.  A  more  versatile  structural  analysis  can  be  used  to  estimate 
the  kriging  equations  and  a  varying  neighborhood  could  be  implemented  in  response  to  range 
fluctuations  or  limit  krigings  use  to  only  isotropic  images.  The  fastest  product  that  would  produce 
high  quality  images  should  use  universal  kriging  with  no  regression  and  incorporate  the  interpolation 
method  for  semi-variogram  estimation. 
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Appendix  A.  Image  Sections 


The  following  is  information  on  the  image  sections  enhanced  for  this  research.  The  data 
includes  the  image  number,  (the  number  used  for  all  following  tables,  semi-variograms,  and  data 
files  output  from  the  program),  information  about  the  image  (e.g.,  urban,  countryside  or  water) 
,  and  the  column  and  row  offsets  input  into  the  program.  An  image  size  of  40  x  40  was  used  in 
all  sub-sampled  images.  In  general  the  SPOT  1  image  consisted  of  countryside  with  a  small  town, 
small  streams,  and  river  valleys.  The  SPOT  2  image  was  of  the  Washington  DC  area  including  the 
Mall  area  and  Potomac  river.  Specifics  of  each  image  are  listed  with  the  following  pixel  data. 


Table  A.l.  SPOT  1,  Channel  2 


SPOT  1 

Description 

Statistics 

Image  Number 

z 

s 

Offset  000,210 

Countryside 

116.19 

40.28 

Offset  108,168 

Countryside 

97.60 

21.74 

Offset  365,370 

Small  Town 

102.52 

46.71 

Offset  390,130 

Countryside 

74.00 

38.65 

Table  A.2.  SPOT  1,  Channel  3 


SPOT  1 

Description 

Statistics 

Image  Number 

z 

5 

Offset  000,210 

Countryside 

102.22 

38.78 

Offset  108,168 

Countryside 

78.38 

21.21 

Offset  365,370 

Small  Town 

90.94 

44.79 

Offset  390,130 

Countryside 

62.51 

34.06 
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Table  A.3.  SPOT  2,  Channel  2 


SPOT  2 

Description 

Statistics 

Image  Number 

z 

s 

Offset  090,250 

Urban 

61.88 

59.83 

Offset  100,130 

Urban 

89.03 

Offset  100,470 

Urban 

124.64 

67.87 

Offset  100,570 

Urban 

139.42 

53.20 

Offset  150,100 

Water 

23.41 

14.77 

Urban 

121.85 

54.76 

Offset  258,000 

Urban 

61.88 

42.74 

Table  A. 4.  SPOT  2,  Channel  3 


SPOT  2 

Description 

Statistics 

Image  Number 

z 

s 

Offset  090,250 

Urban 

67.35 

61.37 

Offset  100,130 

Urban 

99.22 

Offset  100,470 

Urban 

127.33 

68.32 

Offset  100,570 

Urban 

139.12 

55.98 

Offset  150,100 

Water 

83.32 

18.80 

Offset  290.050 

Urban 

130.17 

55.64 

Offset  258,000 

Urban 

28.48 

24.45 
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Appendix  B.  Program  Operation 


The  test  program  was  written  in  Turbo  Pascal  Version  6.0  to  perform  an  einalytical  comparison 
between  kriging  and  cubic  convolution  (current  state  of  the  art).  It  allows  variations  in  trend 
analysis,  model  selection,  neighborhood  size,  sub-sample  size,  and  kriging  method.  Each  selection 
is  performed  by  setting  a  flag  in  the  program. 

B.  1  Flags 

Flags  are  used  within  the  program  to  vary  the  test  conditions. 

B.1.1  Regression.  To  illiminate  the  regression,  the  following  flag  C2in  be  set: 
regression  =  XXXXX; 

where  XXXXX  will  read  “true”  if  regression  will  be  used  and  “false”  if  regression  will  not  be  used. 
For  variations  in  the  regression  terms  within  the  trend  analysis,  the  following  flag  is  set: 

second. order  =  XXXXX; 

where  XXXXX  will  read  “true”  if  second-order  regression  will  be  used  and  “false”  if  second-order 
regression  will  not  be  used.  Second-order  includes  the  eind  terms  into  the  regression  equation 
to  fit  the  polynomial  to  the  data  rather  than  a  linear  function. 

B.1.2  Models.  To  vary  the  model  to  be  used  for  the  semi-variogram,  the  following  flag  is 
set: 

model  =  X; 

where  X  is  set  to  “0”  for  no  model  (interpolation  between  data  points  is  used  for  the  kriging 
equations),  set  to  “1”  to  use  the  spherical  model,  set  to  “2”  for  the  square  root  model,  and  set  to 
“3”  for  the  exponential  model. 

B.1.3  Nugget  Effect.  To  implement  the  nugget  effect,  the  following  flag  is  set: 
nugget  =  XXXXX; 

where  XXXXX  will  read  “true”  if  the  nugget  effect  is  to  be  used  and  “false”  if  is  not  to  be  used. 
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B.1.4  Kriging.  To  vary  the  kriging  method  used  within  the  prograun,  the  following  flag  is 
set: 

kriging  =  1; 

where  X  represents  the  type  of  kriging  used.  The  flag  is  set  to  “0”  if  ordinary  kriging  is  used  and 
“1”  if  universal  kriging  is  used. 


B.2  Image  Data 

The  program  prompts  the  user  to  input  the  name  of  the  image  file  used  for  a  particular 
test.  The  program  lists  all  available  image  files  within  the  directory  specified  in  the  program.  The 
directory  specified  within  the  program  may  be  altered  to  adlow  the  program  to  search  any  directory 
to  accommodate  any  specific  directory  structure. 

After  the  image  file  is  selected,  the  program  prompts  for  a  size  in  width  and  height.  At  this 
point,  the  size  of  the  image  (width  and  height  in  pixels)  section  to  be  studied  will  be  input  at  each 
prompt.  Size  was  determined  from  viewing  the  image  and  selecting  a  proper  section.  Sections  of 
more  than  400  pixels  are  prohibited  due  to  image  file  limitations. 

The  program  also  prompts  for  column  and  row  offsets.  This  is  the  pixel  point  where  the 
upper-right-hand  corner  of  the  selected  section  will  start.  Starting  position,  like  the  size,  was 
determined  from  viewing  the  image  and  selecting  the  starting  position  of  the  proper  section.  The 
upper  right  hand  corner  of  the  original  image  is  the  “0,0”  point. 

The  user  must  be  careful  not  to  position  the  test  section  off  the  original  image  or  include 
sections  of  telemetry  and  calibration  data.  Then,  the  image  section  is  sub-sampled  to  generate  the 
test  file  to  be  kriged.  The  original  image  is  saved  to  a  file  with  an  .ORI  extension.  This  file  will  be 
used  later  to  compare  with  the  enhanced  images. 

B.3  Regression 

The  regression  is  performed  using  procedures  included  in  the  Quinn-Curtis  Huge  Virtual  Array 
Toolbox.  Procedure  HMttliipleReg  uses  least  squares  method  to  choose  the  best-fitting  model  which 
minimizes  the  sum  of  squares  of  the  distances  between  the  observed  responses  and  predicted  by  the 
fitted  model  (22:196).  The  regression  algorithm  takes  the  following  matrix  form  (22:196): 

X  -  matrix  of  independent  data  values, 
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Y  -  vector  of  dependent  variable  values, 

B  -  vector  of  regression  coefficients  to  solve  for. 


XB  =  Y 
X'^XB  =  X'^Y 

B  =  {X'^X)-^Y 

The  vector  B  contains  coefficients  of  the  regression  fit  of  the  data.  McGee  states; 


“An  advantage  of  the  least-squares  polynomial  regression  technique  is  that  it  is  capable 
of  very  closely  modeling  the  trend  . .  .the  trend  is  retained  by  the  polynomial,  and  can 
be  recombined  with  the  kriged  residuals.”  (19:3-7) 

Procedure  HStat Analysis  performs  statistical  analysis  of  the  regression  results  ^lnd  outputs  the 
residual  vector  (22:204).  Regression  output  consists  of  mean,  variance,  standard  deviation  of  the 
original  data  as  well  as  the  resulting  residuals.  The  mean  of  the  residuals  should  be  approximately 
zero  as  a  quick  verification  of  the  calculation.  The  semi-variogram  procedure  uses  the  residual 
vector  supplied  by  the  regression  to  performed  the  necessary  calculations. 

B.4  Semi-Variogram 

Zero  and  ninety  degree  semi-variograms  are  calculated  for  use  in  developing  a  model  for  the 
kriging  equations.  The  semi-variogram  procedure  takes  advamtage  of  the  regularly  grid  structure 
of  the  digital  image  to  simplify  the  calculations. 

B.4  I  Calculations  and  Model  Formulation.  The  variogram  describes  the  variance  of  the 
difference  of  samples  within  the  image  and  is  calculated  as; 

Each  pixel  is  considered  to  be  one  lag  (h)  apart  with  calculations  performed  for  half  of  the  possible 
lags  contained  within  the  image  section  studied. 
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Once  the  semi-variogram  is  calculated,  the  sill  (the  ordinary  variance  of  the  image  data)  is 
calculated  by  the  following: 


«=i 


where 


B.4  S  Data  and  Phi  Generation.  GNUPLOT  macro  files  are  produced  from  the  experimen¬ 
tal  zero  and  ninety  degree  semi-variograms  and  combined  with  a  plot  of  the  sill  and  the  model  if 
desired.  The  following  is  a  typical  macro  file: 


set  term  eepic 

set  output  "kl00-470alcl.tex" 

set  size  1,1 

set  nokey 

set  zrange  [0:10] 

set  yrange  [0:6600] 

set  tics  out 

sot  xlabel  "h'' 

set  ylabel  "$\gainaa(h)$" 

set  title  "Variogram  (Olfset  100,470;  First  Order)" 

plot  3860,  "klOO-470. vOO"  s  linesp,  296.447*8qrt(x)  w  line 


When  run  with  GNUPLOT,  this  macro  will  produce  the  DTgX  file  in  Figure  B.l. 

The  semi-variogram  data  and  plots  were  used  to  determine  the  most  accurate  model  to  use, 
the  level  of  anisotropy  within  the  data,  and  any  variations  in  the  range  of  influence  between  various 
image  sections.  Any  simplifying  assumption  that  can  be  made  with  respect  to  the  semi-variograms 
will  shorten  processing  time  and  decrease  programming  complexity  when  this  is  included  within 
the  image  processing  software  later  (not  performed  for  this  research). 

8.4.3  Output  Files.  The  program  produces  three  files: 

1.  File  one  contains  zero  degree  semi-variogram  data  labeled  as  follows: 

kXXX-YYY.vOO 


2.  File  two  contains  ninety  degree  semi-vanogram  data  labeled  as  follows; 


kXXX-YYY.v90 
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Variogram  (Offset  100,470;  First  Order) 


Figure  B.l.  GNUPLOT  Example. 

3.  File  three  contains  the  GNUPLOT  macro  file  labeled  as  follows; 

kXXX-YYY.gnu 

XXX  represents  the  row  and  YYY  the  column  offset  value. 

B.5  Kriging  Equations 

Once  the  semi-variogram  calculations  are  completed,  both  ordinary  and  universal  kriging 
matrices  were  generated  using  a  16-pixel  neighborhood.  Both  methods  are  used  to  verify  any 
differences  between  the  two.  Ordinary  kriging  is  simpler  to  implement  and  is  less  processor  intensive 
but  universal  kriging  can  handle  any  remaining  local  trend  within  the  data. 

B.  5. 1  Ordinary  Kriging.  The  ordinary  kriging  matrix  is  generated  using  the  chosen  model 
to  estimate  the  semi-variogram  values  for  lag  (/»)  between  the  pixels  within  the  neighborhood.  The 
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ordinary  kriging  matrix  takes  the  following  form: 


1 

1 

1 

7(An) 

7(Ai2)  • 

■  7(Aln) 

7(^21) 

7(^22)  • 

•  7(A2n) 

7(/»nl) 

7(An2)  • 

7(hnn) 

where  n  =  16. 

To  calculate  hij ,  the  distance  between  i  and  j  are  computed  via  the  following  algorithm; 


h  :=  Sqrt(Sqr((i  mod  neighborhood)  -  (j  mod  neighborhood)) 

Sqr((i  div  neighborhood)  -  (j  div  neighborhood))) 

The  value  for  h  is  input  into  the  model  or  interpolation  and  the  result  is  y(hij). 

B.5.2  Universal  Kriging.  The  universal  kriging  matrix  is  generated  like  that  of  the  ordinary 
matrix  with  the  exception  of  the  inclusion  of  term  representing  the  locations  of  the  pixels  within 
the  neighborhood.  To  simplify  the  calculations  and  reduce  the  number  of  matrices,  pixel  one  in  the 
upper  left  hand  corner  of  the  neighborhood  represents  the  "0,0”  pixel  or  the  control  point  of  the 
neighborhood  (Figure  3.1).  The  universal  kriging  matrix  takes  the  following  form; 


0 

0 

0 

1 

1 

1 

0 

0 

0 

Xi 

X2 

Xn 

0 

0 

0 

Vi 

Y2 

Yn 

1 

Xi 

y{hn) 

7(^12)  • 

•  7(^1.) 

1 

X2 

n 

7(h2i) 

7(^22)  • 

•  7(*2r;) 

1 

X„ 

Tn 

7(*».i) 

7(/»n2)  • 

7(^r»n) 

B.5.3  Right-HandSide.  Three  right-hand-side  (RHS)  terms  were  developed  to  correspond 
to  the  three  pixel  estimates  required  within  the  neighborhood.  These  RHS  terms  are  designated 
TOP,  LEFT,  and  MIDDLE  and  Figure  3.1  graphically  displays  the  position  of  pixels  used  to 
create  the  RHS  estimates.  The  RHS  terms  use  the  selected  model  to  estimate  the  semi-variogram 
value  of  the  lag  from  the  estimated  point  to  all  points  in  the  neighborhood  (7(ft„p)).  A  16-pixel 
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neighborhood  will  result  in  a  17-term  RHS  vectors  for  ordin£iry  kriging  and  19-term  RHS  vectors 
for  universal  kriging.  The  following  demonstrate  the  format  for  the  RHS  matrices: 


•  Ordinary  kriging  RHS 


•  Universal  kriging  RHS 


■ 

1 

y(hip) 

7(*2p) 

7(hnp) 


1 

Xp 

Yp 

y{hip) 

y(h2p) 


[  7(hfip)  J 

The  three  RHS  terms  are  combined  with  their  corresponding  kriging  matrices  to  make  a  system  of 
equations  in  matrix  form  to  solve  for  the  weights  in  the  estimator. 


B.5.4  Solution.  The  matrices  and  vectors  are  feed  to  the  Quinn-Curtis  procedure  H  Gauss  Jor¬ 
dan  in  the  matrix  form: 


AX  =  B 


where  A  is  the  square  kriging  matrix,  X  is  the  vector  of  weights,  and  B  is  the  RHS  vector.  A 
matrix  inversion  of  the  matrix  system  takes  place  to  solve  for  the  weights. 

Three  inversions  take  place  that  solve  for  the  weights  at  the  estimate  pixel  locations  and 
these  weights  with  the  neighborhood  are  moved  throughout  the  image  section  to  generate  the 
kriged  image. 

B.6  Cubic  Convolution 

Cubic  convolution  is  implemented  as  discussed  in  Chapter  II  for  a  direct  comparison  with  the 
kriging  techniques. 
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B.  7  Recombination 


The  final  step  of  the  kriging  process  is  the  recombination  of  the  polynomial  equation  with 
the  kriged  residuals.  Both  ordinary  and  universal  kriging  are  compared  to  the  original  image  to 
determine  any  difference  in  techniques. 

B.8  Image  Comparison 

The  image  files  generated  previously  are  use  in  an  analytical  comparison  between  kriging  and 
cubic  convolution  techniques.  The  files  for  kriging  and  cubic  convolution  are  subtracted  for  the 
original  file  to  generate  subtracted  images.  These  may  be  displayed  as  is  but  this  was  not  performed 
for  this  research.  The  difference  of  all  pixels  is  summed  and  a  mean  difference  is  calculated.  A 
similar  calculation  is  performed  to  compute  the  standard  deviation  of  the  data.  Both  the  mean 
and  standard  deviation  for  kriging  and  cubic  convolution  are  written  to  the  display  for  collection. 
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Appendix  C.  Semi-Variograms  of  Residuals 


The  following  are  semi-variograms  generated  from  image  sections  selected.  The  diamonds 
represent  the  0®  direction  calculations  and  the  plus  signs  represent  the  90®  direction  calculations. 
This  data  is  used  to  determine  if  the  range  of  influence  and  anisotropy  vary  between  differing  image 
sections.  Image  sections  from  multiple  channels  of  SPOT  data  were  used  to  see  if  varying  the 
observed  wavelength  causes  any  appreciable  difference  in  the  semi-variogram  calculations.  Similar 
sections  of  differing  channels  are  included  on  the  same  pages  for  ease  of  comparison. 

C.  I  Variograma 


Variogram  (Offset  000,210;  First  Order) 


Figure  CM.  SPOT  1  Channel  2. 


C-1 


Figure  C.2.  SPOT  1  Channel  2. 


Figure  C.3.  SPOT  1  Channel  2. 
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Figure  C.4.  SPOT  1  Channel  2. 


Variogram  (Offset  090,250;  First 
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Figure  C.5.  SPOT  2  Channel  2. 


Figure  C.6.  SPOT  2  Channel  2. 


Figure  C.7.  SPOT  2  Channel  2. 
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Figure  C.8.  SPOT  2  Channel  2. 


Variogram  (Offset  150,100;  First  Order) 
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Figure  C.9.  SPOT  2  Channel  2. 


Figure  C.IO.  SPOT  2  Channel  2. 


Figure  C.H.  SPOT  2  Channel  2. 
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Figure  C.12.  SPOT  2  Channel  3. 


Variogram  (Offset  000,210;  Firs 
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Figure  C.13.  SPOT  1  Channel  3. 


Figure  C.14.  SPOT  1  Channel  3. 


Figure  C.15.  SPOT  1  Channel  3. 
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Figure  C.18.  SPOT  2  Channel  3. 


Figure  C.19.  SPOT  2  Channel  3. 
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C.2  Exponential  Model 


Veiriogram  (Offset  000,210;  First  Order) 


Figure  C.24.  SPOT  1  Channel  2. 


Variogram  (Offset  108,168;  First  Order) 


Figure  C.25.  SPOT  1  Channel  2. 
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Figure  C.26.  SPOT  1  Channel  2. 


Variogram  (Offset  390,130;  Firs 
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Figure  C.27.  SPOT  1  Channel  2. 


Figure  C.31.  SPOT  2  Channel  2. 
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Figure  C.33,  SPOT  2  Channel  2. 
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Figure  C.34.  SPOT  2  Channel  2. 
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Figure  C.36.  SPOT  1  Channel  3. 
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Figure  C.37.  SPOT  1  Channel  3. 


Figure  C.38.  SPOT  1  Channel  3. 
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Figure  C.39.  SPOT  1  Channel  3. 


Figure  C.40.  SPOT  2  Channel  3. 
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Figure  C.41.  SPOT  2  Channel  3. 


Figure  C.42.  SPOT  2  Channel  3. 


C-22 


Figure  C.43.  SPOT  2  Channel  3. 


Figure  C.44.  SPOT  2  Channel  3. 
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Figure  C.45.  SPOT  2  Channel  3. 


Figure  C.46.  SPOT  2  Channel  3. 
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C.3  Spherical  Model 
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Figure  C.50.  SPOT  1  Channel  2. 
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Figure  C.55.  SPOT  2  Channel  2. 


7(/i) 


Figure  C.56.  SPOT  2  Channel  2. 
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Figure  C.57.  SPOT  2  Channel  2. 
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Figure  C.58.  SPOT  2  Channel  3. 


Figure  C.59.  SPOT  1  Channel  3. 
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Figure  C.66.  SPOT  2  Channel  3. 


Figure  C.67.  SPOT  2  Channel  3. 
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Figure  C.68.  SPOT  2  Channel  3. 
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Variogram  (Offset  290,050;  First  Order) 
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Figure  C.69.  SPOT  2  Channel  3. 
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C.4  Nugget  (Exponential  Model). 


Figure  C.71.  SPOT  1  Channel  2. 
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Figure  C.74.  SPOT  2  Channel  2. 


Figure  C.75.  SPOT  2  Channel  2. 
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C-40 


C-41 


Figure  C.80.  SPOT  2  Channel  2. 
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C.5  Nugget  (Spherical  Model). 


Figure  C.81.  SPOT  1  Channel  2. 


Figure  C.82.  SPOT  1  Channel  2. 
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Figure  C.83.  SPOT  1  Channel  2. 
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Figure  C.85.  SPOT  2  Channel  2. 


Figure  C.86.  SPOT  2  Channel  2. 
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Figure  C.88.  SPOT  2  Channel  2. 
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Figure  C.89.  SPOT  2  Channel  2. 
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Figure  C.90.  SPOT  2  Channel  2. 
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Appendix  D.  Test  Statistics  (Mean  and  Standard  Deviation) 


The  following  is  data  generated  in  comparison  of  kriging  and  cubic  convolution  as  well  as 
the  spherical  and  exponential  models  with  interpolation  method.  Images  were  sub-sampled  and 
data  regenerated  using  all  techniques.  The  reproduced  image  was  then  subtracted  from  the  original 
data,  a  sample  mean  (i)  and  standard  deviation  (s)  of  the  difference  was  calculated  for  all  pixels, 
and  this  data  was  compared  between  the  two  techniques. 

D.l  Kriging  and  Cubic  Convolution  Comparison 

If  the  tests  statistics  indicate  that  Kriging  is  the  better  of  the  two  means  or  variances,  a  K 
will  be  entered  in  the  table  under  the  column  Best,  a  C  will  be  entered  if  Cubic  convolution  is  best 
and  a  -  if  neither  is  considered  bets  by  the  statistics. 


Table  D.l.  Comparison  Data  of  Kriging  and  Cubic  Convolution,  Channel  2. 


SPOT  1 

Kriging 

Cubic  Convolution 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

K 

X 

S 

Offset  000,210 

16.34 

17.27 

16.64 

17.94 

.35 

.93 

Offset  108,168 

ItBHI 

8.77 

12.35 

9.78 

2.68 

.80 

la 

Offset  365,370 

16.63 

■mi 

.34 

.67 

- 

■a 

Offset  390,130 

12.28 

11.04 

13.86 

12.25 

2.82 

.81 

o 

B 

Table  D.2.  Comparison  Data  of  Kriging  and  Cubic  Convolution,  Channel  3. 


SPOT  1 

Kriging 

Cubic  Convolution 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

B 

Bl 

Offset  000,210 

19.09 

19.31 

19.43 

19.9 

.36 

.93 

- 

- 

Offset  108,168 

12.75 

9.59 

14.3 

10.36 

3.23 

.85 

B 

Bl 

Offset  365,370 

22.35 

18.08 

22.5 

18.32 

.17 

.97 

- 

- 

Offset  390,130 

14.96 

12.38 

16.43 

13.9 

2.3 

.78 

B 

B 
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Table  D.3.  Comparison  Data  of  Kriging  and  Cubic  Convolution,  Channel  2. 


SPOT  2 

Kriging 

Cubic  Convolution 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

5 

t 

F, 

El 

m 

Offset  090,250 

18.95 

19.25 

19.56 

19.14 

.66 

1.01 

- 

- 

Offset  100,130 

22.65 

19.79 

23.19 

20.25 

.56 

.95 

- 

- 

Offset  100,470 

22.67 

21.35 

24.12 

21.81 

1.39 

.95 

- 

- 

Offset  100,570 

34.01 

27.96 

35.83 

30.75 

1.28 

.82 

- 

Q 

Offset  150,100 

12.32 

9.56 

13.68 

10.66 

2.79 

.80 

Q 

Q 

Offset  290,050 

27.37 

26.66 

28.27 

28.65 

.67 

.87 

- 

- 

Offset  258,000 

15.98 

12.94 

16.5 

13.9 

.81 

.87 

- 

- 

Table  D.4.  Comparison  Data  of  Kriging  and  Cubic  Convolution,  Channel  3. 


SPOT  2 

Kriging 

Cubic  Convolution 

Test  Statistics 

Best 

Image  Number 

X 

8 

X 

S 

t 

Fo 

a 

5 

Offset  090,250 

20.62 

20.3 

21.45 

20.34 

.85 

.99 

- 

- 

Offset  100,130 

25.13 

21.36 

25.91 

21.51 

.98 

- 

- 

Offset  100,470 

24.19 

22.84 

25.61 

23.33 

1.28 

.95 

- 

- 

Offset  100,570 

36.62 

29.3 

38.66 

32.13 

1.38 

.83 

- 

Q 

Offset  150,100 

16.11 

14.19 

18.12 

15.95 

2.77 

.79 

o 

O 

Offset  290,050 

28.07 

27.66 

28.93 

30.03 

.62 

.84 

- 

m 

Offset  258,000 

17.2 

14.8 

T:inil 

15.96 

1.33 

.85 

- 

- 
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D.2  Interpolation  and  Spherical  Model  Comparison 


If  the  tests  statistics  indicate  that  interpolation  is  the  better  of  the  two  means  or  variances, 
an  1  will  be  entered  in  the  table  under  the  column  Best,  a  S  will  be  entered  if  the  Spherical  model 
is  best  and  a  -  if  neither  is  considered  bets  by  the  statistics. 


Table  D.5.  Comparison  Data  of  Interpolation  and  Spherical  Model,  Channel  2. 


SPOT  1 

Spherical  Model 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

F„ 

B 

m 

Offset  000,210 

16.32 

17.25 

16.34 

17.27 

.02 

.99 

- 

- 

Offset  108,168 

11.08 

8.72 

11.15 

8.77 

.16 

.98 

- 

- 

Offset  365,370 

20.38 

16.67 

20.22 

16.63 

-.20 

1.0 

- 

- 

Offset  390,130 

12.75 

10.92 

12.28 

11.04 

-.89 

.97 

- 

- 

Table  D.6.  Comparison  Data  of  Interpolation  and  Spherical  Model,  Channel  3. 


SPOT  1 

Spherical  Model 

Interpolation 

Test  Statistics 

Image  Number 

X 

s 

X 

s 

t 

Fo 

B 

m 

Offset  000,210 

19.23 

19.21 

19.09 

19.31 

-.28 

.95 

- 

- 

Offset  108,168 

12.67 

9.59 

.17 

.99 

- 

- 

Offset  365,370 

22.56 

17.9 

22.35 

18.08 

-.24 

.98 

- 

- 

Offset  390,130 

14.98 

12.27 

14.96 

12.38 

-.25 

.98 

- 

- 

Table  D.7.  Comparison  Data  of  Interpolation  and  Spherical  Model,  Channel  2. 


SPOT  2 

Spherical  Model 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

B 

m 

Offset  090,250 

19.69 

19.01 

18.95 

19.25 

-.81 

.97 

- 

- 

Offset  100,130 

23.02 

19.66 

22.65 

19.79 

-.39 

.98 

- 

- 

Offset  100,470 

23.32 

21.31 

22.67 

21.35 

.99 

- 

- 

Offset  100,570 

35.03 

28.03 

34.01 

27.96 

-.75 

1.0 

- 

- 

Offset  150,100 

12.22 

9.23 

12.32 

9.56 

.22 

.93 

- 

- 

Offset  290,050 

28.36 

26.55 

27.37 

26.66 

-.77 

.99 

- 

- 

Offset  258,000 

16.19 

12.9 

15.98 

12.94 

.-.34 

.99 

- 

- 

Table  D.8.  Comparison  Data  of  Interpolation  and  Spherical  Model,  Channel  3. 


SPOT  2 

Interpolation 

Spherical  Model 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

8 

t 

Fo 

Q 

m 

Offset  090,250 

21.27 

20.09 

20.62 

20.3 

-.47 

.92 

- 

- 

Offset  100,130 

25.5 

EBI 

25.13 

-.36 

.98 

- 

- 

Offset  100,470 

24.91 

1^3 

22.84 

-.65 

.99 

- 

- 

Offset  100,570 

37.76 

29.23 

36.62 

29.3 

-.81 

.99 

- 

- 

Offset  150,100 

15.92 

14.16 

16.11 

14.19 

.27 

.99 

- 

- 

Offset  290,050 

28.86 

27.57 

28.07 

27.66 

-.59 

.99 

- 

- 

Offset  258,000 

17.43 

14.74 

17.2 

14.8 

-.32 

.99 

-  j 

- 
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D.3  Interpolation  and  Exponential  Model  Comparison 


If  the  tests  statistics  indicate  that  interpolation  is  the  better  of  the  two  means  or  variances, 
an  I  will  be  entered  in  the  table  under  the  column  Best,  an  E  will  be  entered  if  the  Exponential 
model  is  best  and  a  -  if  neither  is  considered  bets  by  the  statistics. 


Table  D.9.  Comparison  Data  of  Interpolation  and  Exponential  Model,  Channel  2. 


SPOT  1 

Exponential  Mode 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

B 

m 

Offset  000,210 

16.34 

17.26 

16.34 

17.27 

0.0 

.99 

- 

- 

Offset  108,168 

11.06 

8.69 

11.15 

8.77 

.21 

.98 

- 

- 

Offset  365,370 

20.42 

16.72 

20.22 

16.63 

-.24 

1.0 

- 

- 

Offset  390,130 

12.74 

10.89 

12.28 

11.04 

-.87 

.97 

- 

- 

Table  D.IO.  Comparison  Data  of  Interpolation  and  Exponential  Model,  Channel  3. 


SPOT  1 

Exponential  Model 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

B 

m 

Offset  000,210 

19.25 

19.21 

19.09 

19.31 

-.17 

.98 

- 

- 

Offset  108,168 

12.63 

9.55 

12.75 

9.59 

.26 

.99 

- 

- 

Offset  365,370 

22.6 

18 

22.35 

18.08 

-.28 

.99 

- 

- 

Offset  390,130 

14.96 

12.24 

14.96 

12.38 

0.0 

.98 

- 

- 
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Table  D.ll.  Comparison  Data  of  Interpolation  eind  Exponential  Model,  Channel  2. 


SPOT  2 

Exponential  Model 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

5 

t 

Fo 

B 

m 

Offset  090,250 

19.76 

19.09 

18.95 

19.25 

-.87 

.98 

- 

- 

Offset  100,130 

23.11 

19.66 

22.65 

19.79 

-.48 

.98 

- 

- 

Offset  100,470 

23.4 

21.35 

22.67 

21.35 

-.71 

1.0 

- 

- 

Offset  100,570 

35.04 

27.99 

34.01 

27.96 

-.77 

1.0 

- 

- 

Offset  150,100 

12.19 

9.19 

12.32 

9.56 

.28 

.92 

- 

- 

Offset  290,050 

28  43 

26.51 

27.37 

26.66 

-.83 

.98 

- 

- 

Offset  258,000 

16.23 

12.91 

15.98 

12.94 

-.40 

.99 

- 

- 

Table  D.12.  Comparison  Data  of  Interpolation  and  Exponential  Model,  Channel  3. 


SPOT  2 

Exponential  Model 

Interpolation 

Test  Statistics 

Best 

Image  Number 

X 

8 

X 

8 

t 

Fo 

B 

m 

Offset  090,250 

21.36 

20.12 

20.62 

20.3 

-.76 

98 

- 

- 

Offset  100,130 

25.57 

21.25 

25.13 

21.36 

-.43 

.98 

- 

- 

Offset  100,470 

25 

22.78 

24.19 

22.84 

-.73 

.99 

- 

- 

Offset  100,570 

37.78 

29.19 

36.62 

-.82 

.99 

- 

- 

Offset  150,100 

15.86 

14.18 

■CTfl 

14.19 

.36 

.99 

- 

- 

Offset  290,050 

28.95 

27.48 

28.07 

27.66 

-.66 

.99 

- 

- 

Offset  258,000 

17.48 

14.73 

17.2 

14.8 

-.39 

.99  1  - 

- 
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D.4  Non-Nugget  and  Nugget  Comparison 


0.4.1  Spherical  Model.  If  the  tests  statistics  indicate  that  non-nugget  is  the  better  of  the 
two  means  or  variances,  an  NN  will  be  entered  in  the  table  under  the  column  Best,  a  N  will  be 
entered  if  the  Exponential  model  is  best  and  a  -  if  neither  is  considered  bets  by  the  statistics. 


Table  D.13.  Comparison  Data  of  Non-Nugget  and  Nugget,  Channel  2,  Spherical  Model. 


SPOT  1 

Non-Nugget 

Nugget 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

X 

m 

Offset  000,210 

16.32 

17.25 

16.83 

17.64 

.60 

.95 

- 

- 

Offset  108,168 

11.08 

8.72 

11.04 

8.58 

-.09 

1.0 

- 

- 

Offset  365,370 

20.38 

16.67 

22.12 

2.09 

.87 

NN 

- 

Offset  390,130 

12.75 

10.92 

13.11 

10.75 

.69 

1.0 

- 

- 

Table  D.14.  Comparison  Data  of  Non-Nugget  and  Nugget,  Channel  2,  Spherical  Model. 


SPOT  2 

Nugget 

Non-Nugget 

Test  Statistics 

Best 

Image  Number 

X 

« 

X 

s 

t 

Fo 

X 

s 

Offset  090,250 

19.01 

20.45 

3.06 

.86 

NN 

- 

Offset  100,130 

19.66 

24.46 

1.51 

.97 

- 

- 

Offset  100,470 

21.31 

25.69 

2.27 

.93 

NN 

- 

Offset  100,570 

35.03 

28.03 

39.15 

29.44 

2.99 

.90 

NN 

- 

Offset  150,100 

12.22 

9.23 

11.65 

8.19 

-1.36 

1.27 

- 

m 

Offset  290,050 

28.36 

26.55 

31.5 

26.96 

2.44 

.96 

NN 

- 

Offset  258,000 

16.19 

12.9 

17.63 

13.54 

2.27 

.91 

- 

- 
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D.5  Non- Nugget  and  Nugget  Comparison 


Table  D.15.  Comparison  Data  of  Non- Nugget  and  Nugget,  Channel  2,  Exponential  Model. 


SPOT  1 

Non-Nugget 

Nugget 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

s 

t 

Fo 

X 

s 

Offset  000,210 

16.34 

17.26 

16.9 

17.7 

.66 

.95 

- 

- 

Offset  108,168 

11.06 

8.69 

11.05 

8.5 

-.02 

1.03 

- 

- 

Offset  365,370 

20.42 

16.72 

22.28 

18 

2.22 

.86 

NN 

- 

Offset  390,130 

12.74 

10.89 

13.14 

10.73 

.77 

1.03 

- 

- 

Table  D.16.  Comparison  Data  of  Non-Nugget  and  Nugget,  Channel  2,  Exponential  Model. 


SPOT  2 

Non-Nugget 

Nugget 

Test  Statistics 

Best 

Image  Number 

X 

s 

X 

5 

t 

Fo 

s 

Offset  090,250 

19.76 

19.09 

22.97 

20.54 

3.37 

.86 

NN 

- 

Offset  100,130 

23.11 

19.66 

24.88 

19.98 

1.85 

.96 

- 

- 

Offset  100,470 

23.4 

21.35 

25.92 

22.07 

2.42 

.93 

NN 

- 

Offset  100,570 

35.04 

27.99 

39.23 

29.53 

3.03 

.89 

NN 

- 

Offset  150,100 

12.19 

9.19 

11.63 

8.16 

-1.34 

1.26 

- 

EM 

Offset  290,050 

28.43 

26.51 

31.81 

27.03 

.96 

NN 

- 

Offset  258,000 

16.23 

12.91 

17.82 

13.65 

BEI 

.89 

NN 

- 
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D.5.1  Exponential  Model.  Appendix  E.  Program 


E.l  Program  Listing 

The  following  is  the  Pascal  code  used  in  thesis  research. 


Uses  CRT.Misc,TP030, 
Support , 
HinHax, 
hmulreg, 
haataath , 
■atasc , 
stats , 
leq. 

Printer; 


type 

line  =  array  [0..799]  ol  byte; 
const 

regression  ~  false; 
second. order  =  false; 

■odel  »  0;  {0»Hone,l“Spherical,2»Square  Root,3’*E3q>onential} 

nugget  B  false; 

neighborhood  «  4; 

kriging  -  1;  {0«>Siaple,l=Universal} 
sub_sa]q)le  °  2 ; 
latex  ■  true ; 

debug  false; 


Tar 

pixel_c,pixel_k 

pixel.o, 

nuB.iv , din , start 
result ,rl,r2, 
rdia.cdia 
roes, cols, 
h.i.j 

■ax.lag, range 

col.of f set ,ro¥_of f set , 

k,kh, pixels, n 

■ean_r es , std_res , Tar.r es , 

■ean,suaz, 

std.deT , variance , 

gc,nug, 

d , det , cond, aaxerr , 

dif_c,dif_k, 

biae_c,bias_k, 

suB_c ,  8ua_c2 ,  stin.k ,  8ua_k2 , 

■ean_c .aeem.k, 

stdT.c , stdT.k , 

z.value 

cpix,kpix, 

reg 


:  double; 

:  byte; 

:  integer ; 

:  longint ; 
:  word; 

;  longint ; 


:  double; 
:  double; 
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gauaO  ,gana90 
cn,kn, 

■M.sT.vL 
dist, 
pixel , 

data,coef .zest .resid, 

gaama , gcoei .res , 

AMat.ANatInv. 

TRHS.NRHS.LRHS. 

TSol.MSol.LSol. 

KSol.CSol 

rst 

fni.fno 

buffer 

fi 

fo 

fbo.fbc.fbk 

si.sj.di.dj.mi.aj 


:  array  [0..4]  of  double; 

:  array  [0 . . 200]  of  double ; 

:  array  [0.  .neighborhood- 1,0. .neighborhood- 1]  of  double; 


:  ARRDESP; 

:  regstats; 

:  string; 

:  line; 

:  file  of  line; 
:  text ; 

:  file  of  byte; 
:  integer; 


Function  Variogran(diat  :  double)  :  double; 
var 

lover .upper  ;  byte; 
begin 

case  nodel  of 
0  :  begin 

if  Trunc(di8t)  =  dist  then 
Variograa  :»  gannaO[Tnmc(dist)] 
else 
begin 

lover  :=  Trunc(dist); 
upper  I*  lover  +  1; 

Variogran  :=  gauaO  [lover] 

+  (gauaO [upper]  -  gannaO [lover] )* (dist  -  lover); 

end;  {else} 
end;  {No  Model} 

1  :  begin 

Variogran  :»  gc*(1.5-0.S*Sqr(dist/range))*(dist/range)  +  nug; 
end;  {Spherical  Model} 

2  :  begin 

Variogran  gc*sqrt (dist)  nug; 

end;  {Square  Root  Model} 

3  :  begin 

Variogran  :■  gc*(l-exp(-dist/range))  +  nug; 
end;  {Exponential  Model} 
end;  {case} 

end;  {Function  Variogran} 


BEGIN 


{  Input  inage  } 

ClrScr ; 

f ni  : ■  Select_File( ' Inage ’ , ' \IMAGES\DATA\* . * ’ , ' SPOTl . CHI’, 1,1, IB, S); 
A88ign(fi, ’\IMAGES\DATA\’+fni) ; 

Reset(f i) ; 

GotoXYd.S); 
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{  Initiali2e  regression  } 

Hrite( ’Inage  width?  ’); 

Read(cols) ; 

Vrite( ’Inage  height?  ’); 

Read(rows) ; 

pixels  rows  *  cols; 
pixel  HVectDef (pixels) : 

Vrite(’Colunn  offset?  ’); 

Read(col_off set) ; 

Write (’Ron  offset?  ’); 

Read(row_off set) ; 

fno  :■  ’K’  +  ThreeDigit(row_offset) 

+  ’-’  +  ThreeDigit(col_offset) ; 

Write  (’Reading  inage _  ’); 

k  0; 

for  i  :»  0  to  rows-l  do 
begin 

Seek(fi,row_offset-»'Sub_aanple*i) ; 

Read(fi, buffer) ; 

Nark.Tine ; 

for  j  :»  0  to  cols-l  do 
begin 

HVectWrtEl (pixel, k,buffer[col_offset-t-sub_sanple*j]  .result) ; 
k  k  +  1; 
end;  {for  j} 
end;  {for  i} 

Writeln; 
if  debug  then 

HNatSaweASC(pixel,fno+’ .pix’ ,9, 3, result) ; 
nean  HVectNean (pixel, v_col,0,sunx) ; 

std.dew  :=  H?ectStDeT(pixel, w.col.O, nean, variance) ; 

Writeln( ’Original  Inage ’ ) ; 

WritelnC  Nean  =  ’,nean:7:2); 

Writeln(’  St  Dev  »  ’ ,std_dev:7;2) ; 

Writeln(’  Variance'  ’ .variance :7: 2) ; 
if  sub.sanple  >  1  then 
begin 

Write ( ’Writing  original  file  before  subsanpling. . . ’) ; 

Seek(f i,ros_offset) ; 

Assign (f bo, fno+’ .ORI’) ; 

Rewrite (f bo) ; 

for  i  0  to  (sub.sanplevrows)-!  do 
begin 
Nark.Tine; 

Read(fi, buffer) ; 

for  j  3  to  (sub_sanple*cols)-l  do 
Writevfbo,buffer[col_offset+j]) ; 
end;  {for  i} 

Close (f bo) ; 
end;  {if} 

Close(f i) ; 
if  regression  then 
begin 

if  second_order  then 
nun_iv  :»  4 
else 
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nuM_iT  2; 

data  HMatDef  (pixels, niui_lT) ; 

zest  HVectDel (pixels) ; 

resid  HVectDef (pixels) ; 

coel  HVectDef  (nuit_iT+l)  ; 

{  Initialize  X  independent  variable  } 
k  0; 

Mark.Tine ; 

for  i  0  to  rows-l  do 

for  j  0  to  cols-1  do 
begin 

HNatVrtEl(data,k,0,j .result) ; 
k  k  +  1; 

end;  {for  j} 

{  Initialize  Y  independent  variable  } 
k  0; 

Mark.Tine ; 

for  i  :■  0  to  rovs-l  do 
for  j  :=  0  to  cols-l  do 
begin 

HMatVrtEl(data,k,l ,i, result) ; 
k  :=  k  +  1; 
end;  {for  j} 
if  second_order  then 
begin 

{  Initialize  X*2  independent  veuriable  } 
k  0; 

Nark.Tine; 

for  i  :=>  0  to  rovs-l  do 
for  j  :=•  0  to  cols-l  do 
begin 

HMatVrtEl(data,k,2,Sqr(j) .result) ; 
k  k  +  1; 
end;  {for  j} 

{  Initialize  Y*2  independent  variable  } 
k  0; 

Hark.Tine ; 

for  i  :*  0  to  rows-l  do 
for  j  :■  0  to  cols-l  do 
begin 

HMatHrtEl(data,k,3,Sqr(i) .result) ; 
k  k  +  1; 
end;  {for  j} 
end;  {if  second. order} 

Vriteln; 

if  second.order  then 

tfrite(’ Calculating  second-order  regression. . . ’) 
else 

Write ( ’Calculating  first-order  regression. . . ’) ; 

HMult ipleReg (data .pixel , coef , nil , TRUE , nil .result ) ; 
for  i  :■  0  to  nun.iv  do 

c{i]  HVectReadEKcoef.i, result)  ; 

Writeln; 

Hrite( ’Calculating  residuals. . . ’) ; 

HStat Analysis (data, pixel, coef .nil .nil, 0.96. nil .resid, zest, rst) ; 
■ean_res  HVectHean(resid,v_col,0,suax) ; 
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std_r«8  HTectSt0ev(reBid,T_col,0,Kean_res,Tar_re8} ; 

Writeln; 

VritelnC ’ Re8idual8 ’ ) ; 

Vriteln(’  Mean  »  ’ ,nean_re8:7:2)  ; 

Uriteln(’  St  Dev  »  ’ ,8td_re8:7:2) ; 

VritelnC’  Variance  •  ’ ,Tar_re8;7:2) ; 
end  {if  regrenaion} 
elae 
begin 
Vriteln; 

end;  {elae  no  regreaaion} 
nax.lag  :>  IMinCcola  div  2, rose  div  2); 
range  :=  5; 

{  Zero  degree  variogran  } 

tfrite( ’Calculating  zero  degree  rariogran. . . ; 
gaanaOCO]  :»  0; 
for  h  :>•  1  to  nu.lag  do 
begin 
Mark.Tine ; 
n  :=  0; 

gannaO [h]  : ~  0 ; 

for  i  :■  0  to  row8-l  do 

for  j  :»  0  to  (cola  -  h  -  1)  do 
begin 

n  n  +  1; 
k  :»  i«cols  +  j; 
kb  :»  i*col8  +  j  +  h; 
if  regreasion  then 
gaanaOCh]  gaMaOCh] 

Sqr(HVectReadEl(rejid,k,rl)  -  HVectRead£l(re8id,kh,r2)) 

elae 

gauaO  [h]  :  >  gauaO  [h] 

-t*  Sqr(HVectReadEl(pizel,k,rl)  -  HVectReadEl(pizel,kh,r2)) ; 

end;  {for  j} 

ganaaO  [h]  ; »  gauaO  [h]  /  (2*n)  ; 
end;  {for  h} 

A88ign(fo,fno+’ .VOO’) ; 

Rewrite (fo) ; 

for  h  :>  0  to  naz.lag  do 

Vriteln(fo,h:3,gania0[h] :  10:2)  ; 

CloaaCfo) ; 

Hriteln; 

{  linety  degree  wariograa  } 

Hrite( ’Calculating  ninety  degree  wariogran. . . ; 
gannaSO [0]  : ■  0 ; 
for  h  1  to  nu.lag  do 
begin 
Kark.Tine ; 
n  0; 

geunia90[b]  :•  0; 

for  j  0  to  cola-l  do 

for  i  :■  0  to  (rowa-h-l)  do 
begin 

n  :■  n  +  1; 
k  :■  iwcola  +  j; 
kb  (i+h)*col8  +  j; 
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if  regression  then 

gannaSO  [h]  :  >  gauaSO  [h] 

+  Sqr(HVectRead£l(resid,k>rl)  -  HVectReadEl(resid,kh,r2)) 

else 

gannaSO  [h]  :  gaiuaRO  [h] 

+  Sqr(HVectReadEl(pixel,k.rl)  -  HVectRead£l(pixel,kh,r2)) ; 

end;  {for  i} 

gannaRO  [h]  :  >  gauaSO  [h]  /  (2*n)  ; 
end;  {for  h} 

A88ign(fo,fno+’ .V90’) ; 

Rewrite (fo) ; 

for  h  :=•  0  to  nu.lag  do 

Vriteln(fo,h:3,ganna90[h3  :10:2); 

Close(fo) ; 

Writeln; 

if  regression  then 
begin 

HArrFree(data) ; 

HArrFreeCzest) ; 

HArrFree(coef) ; 
end;  {if  regression} 

{  Compute  Model  Coefficients  (Zero  Degree)  } 
if  model  >  0  then 
begin 

num_iv  :«  1; 

data  HMatDef  (range-^l.num.iv)  ; 

gamma  HVectDef  (range-*-!)  ; 

gcoef  :=•  HVectDef  (num.iw-t-l)  ; 

{  Initialize  matrices  } 
for  i  :=  0  to  range  do 

HTectWrtEl (gamma, i,g2aima0[i]  result)  ; 
for  i  :>■  0  to  range  do 
case  model  of 

1  :  HVectVrtEl(data,i,(1.6  -  0.5*Sqr(i/range))*(i/range) .result) ; 

2  :  HVectWrtEl(data,i,Sqrt(l) .result) ; 

3  :  HVectHrtEl(data,i,l-Ezp(-i/range) .result) ; 
end;  {case} 

if  nugget  then 
begin 

HMult ipleReg (data , gamma , gcoef , nil , TRUE .nil , result ) ; 
nug  : *■  HVectReadEl (gcoef , 0  .result)  ; 
end  {if} 
else 
begin 

HMultipleReg (data , gamma .gcoef .nil .FALSE .nil .result) ; 
nug  0; 
end;  {else} 

gc  ;•  HVectReadEl (gcoef ,1, result) ; 
end;  {if  model  >  0} 

{  Generate  GMU  Plot  file  } 

A88ign(fo,fno+’ .GHU’) ; 

Rewrite (fo) ; 
if  latex  then 
begin 

Hriteln(fo, 'set  term  eepic’); 

Vriteln(fo, ’set  output  " ’ ,fno, ’2.tox" ’) ; 
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end  {if} 
else 

Vritelndo, ’set  tern  vgalib’); 

Writeln(fo, ’set  size  1,1’); 

Vriteln(fo, ’set  nokey’); 

VritelnCfo, ’set  zrange  [0: ’ ,naz_lag, ’] ’) ; 

Vriteln(fo, ’set  yrange  [0:5500]’); 

WritelnCfo, ’set  tics  out’); 

VritelnCfo, ’set  zlabel  "h"’); 

VritelnCfo, ’set  ylabel  "$\gannaCh)$"’) ; 

Vrite  Cfo,’set  title  "Variogran  COffset  ’, 

ThreeDigit  Cron.of f set) , ’ , ’ , 

ThreeDigitCcol.offset)) ; 
if  regression  then 
if  second.order  then 

VritelnCfo,’;  Second  Order)"’) 
else 

VritelnCfo,’;  First  Order)"’) 

else 

VritelnCfo,’;  Hone)"’); 
if  regression  then 

Vrite  Cfo,’plot  ’  ,»su:_reo:8:2,  ’ ,  "k’, 

ThreeDigit Cros.off set) , , 

ThreeDigit Ccol_off set) ,’ .vOO"  w  linesp,  "k’, 
ThreeDigitCrow.offset) , , 

ThreeDigit Ccol.off set), ’.»90"  a  linesp’) 

else 

Vrite  Cfo,’plot  ’, variance •.8:2,  ’ ,  "k’, 

ThreeDigit Croe.of feet) , , 
ThreeDigitCcol.offset) , ’ .vOO"  w  linesp,  "k’, 
ThreeDigit Croe_off set) , , 

ThreeDigit Ccol_off set) ,’ .v90"  w  linesp’); 
if  nodel  >  0  then 
begin 

VriteCfo,’,  ’,gc:8:3); 
case  model  of 

1  :  VriteCfo, ’eCl. 6*Cx/’, range, ’)-0. 5*Cx/’, range, ’)**3)’) 

2  :  VriteCfo, ’esqrtCx) ’) ; 

3  :  VriteCfo, ’*Cl-oxpC-x/’ , range, ’))’) ; 
end;  {case} 

if  nugget  then 
VritelnCfo,’  +  ’,nug:8:3) 
else 

VritelnCfo) ; 
end  {if} 
else 

VritelnCfo) ; 

CloseCfo) ; 

HArrFreeCdata) ; 

HArrFreeCganna) ; 

HArrFreeCgcoef ) ; 

{  Set  up  Kriging  Equations  } 
if  kriging  ■  0  then 

start  :•■  1  {Simple  kriging} 
else 
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start  :**  3;  {Ibiversal  kriging} 
dia  :»  Sqr (neighborhood)  +  start; 

ANat  HMatDef (diH.din) ; 

TRHS  HVectDel(din); 

MRHS  HVectDef(diii); 

LRHS  HVectDeKdin); 

TSol  :>  HVectDef(din); 

NSol  HVectDef(din): 

LSol  HVectOef(din); 

ANatInv  HMatDef (din,diH) ; 
lor  i  :■  0  to  start-1  do 
lor  j  0  to  start- 1  do 

HNatVrtEl(AMat,i, j ,0, result) ; 
lor  i  :■  start  to  din-1  do 
begin 

HMatHrtEl(AMat,0,i,l, result) ; 

HHatHrtEl ( AMat , i , 0 , 1 .result ) ; 
end;  {lor  i> 
il  kriging  ~  1  then 
begin 
n  ;»  2; 

lor  i  1  to  neighborhood  do 
lor  j  1  to  neighborhood  do 
begin 

n  n  +  1; 

HMatVrtEl (AMat, n,l,i, result) ; 

HMatWrtEl(AHat,n,2,j .result) ; 

HMatHrtEKAHat .  1  .n.i. result) ; 

HHatVrtEK AKat , 2  .n .  j  .result) ; 
end;  {lor  j} 
end;  {else} 

lor  i  :=  0  to  din- (start  +  1)  do 
lor  j  :»  0  to  din- (start  +  1)  do 
begin 

d  Sqrt(Sqr((i  nod  neighborhood)  -  (j  nod  neighborhood)) 

+  Sqr((i  div  neighborhood)  -  (j  div  neighborhood))); 
HMatWrtEl(ANat.i-*-steu:t.  j-*-start.Variogran(d)  .result) ; 
end;  {lor  j} 

{  Set  up  Middle  RHS  lor  kriging  equations  } 

HVectVrtEl (MRHS. 0.1. result) ; 
il  kriging  =  1  then 
begin 

HVectHrtEl (MRHS .1.2. 5 .result) ; 

HVectVrtEl (MRHS .2.2.5. result ) ; 
end; 

k  :»  start; 
lor  i  :■  1  to  4  do 
lor  j  ;•  1  to  4  do 
begin 

d  :■  Sqrt(Sqr(i  -  2.5)  +  Sqr(j  -  2.5)); 

HVectHrtEl (MRHS. k.Variogran(d) .result) ; 

k  k  +  1; 

end; 

HGaussJordan(AMat .MRHS .MSol .AMatlnv .det .cond.naxerr. result) ; 
il  debug  then 
begin 
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HMatSaveASC(ANatI&T,fno+’ .AMI’ ,12, 6, result) ; 
HNatSaTeASC(NRHS,fno+’ .RHM’ , 12, 6, result) ; 
HMatSaTeASC(MSol,fno+>  .HTMM2, 6, result)  ; 
end;  {if  debug} 

{  Set  up  Top  RHS  for  kriging  equations  } 

HVectWrtEKTHHS, 0,1, result) ; 
if  kriging  =  1  then 
begin 

HVectWrtEl (TRHS ,1,2.0, result ) ; 

HVectVrtEl (TRHS ,2,2.5, result ) ; 
end; 

k  start; 
for  i  1  to  4  do 
for  j  :=  1  to  4  do 
begin 

d  :»  Sqrt(Sqr(i  -  2.0)  +  Sqr(j  -  2.5)); 

HVectWrtEl  (TRHS,  k,V2a'iogran(d)  , result)  ; 

k  k  +  1; 

end; 

HGauss Jordsm  (ANat , TRHS ,TSol , AMat Inv , det , cond , nazerr , result ) ; 
if  debug  then 
begin 

HMatSaveASC(TRHS,fno+’  .RHTM2, 6, result)  ; 
HMatSaTeASC(TSol,fno>’ .WTT’ , 12, 6, result) ; 
end;  {if  debug} 

{  Set  up  Left  RHS  for  kriging  equations  } 

HVectWrtEKLRHS, 0,1, result) ; 
if  kriging  =  1  then 
begin 

HVectWrtEl (LRHS ,1,2.5, result ) ; 

HVectWrtEl (LRHS ,2,2.0, result ) ; 
end; 

k  :=  start; 
for  i  :«■  1  to  4  do 
for  j  :=  1  to  4  do 
begin 

d  :■  Sqrt(Sqr(i  -  2.5)  +  Sqr(j  -  2.0)); 

HVectWrtEl (LRHS, k,Variogran(d) , result) ; 

k  k  +  1; 

end; 

HGauss Jordan (ANat ,LRHS ,LSol , ANatInv , det , cond, Baxerr, result) ; 
if  debug  then 
begin 

HNatSaTeASC(LRHS,fno>’ .RHL ’,12, 6, result) ; 
HNatSaveASC(LSol,fno+’ .WTL’ ,12, 6, result) ; 
end;  {if  debug} 

HArrFree(ANat) ; 

HArrFree (ANatInv) ; 

HArrFree(LRHS) ; 

HArrFree (TRHS) ; 

HArrFree (NRHS) ; 
k  0; 

for  i  0  to  neighborhood- 1  do 
for  j  0  to  neighborhood- 1  do 
begin 

«L[i, j]  ;>HVectReadEl(LSol,k+start,rl) ; 
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wT[i , j]  : »HVectReadEl (TSol ,k+8tart ,rl) ; 

«NCi,  j]  :«HVectReadEl(IISol,k-*-8tart,rl) ; 
k  :•  k  +  1; 
end;  {for  j} 

HArrFree(LSol) ; 

HArrFree(TSol) ; 

HAixFree  (NSol) ; 

rdin  :•  roes  *  sub.sample; 

cdin  cols  •  sub.Banple; 

KSOL  HNatDef (rdiB,cdim) ; 

CSOL  HIIatDef  (rdin.cdim)  ; 
for  si  :=  0  to  rdi*-l  do 
for  sj  :>•  0  to  cdin-1  do 
begin 

■i  si  nod  sub_sanple;  di  si  div  sub.Scaple; 

■j  ;»  sj  mod  sub.sanple;  dj  sj  dir  sub.seimple; 

if  (mi  =■  0)  and  (mj  =  0)  then 
{•  Read  original  image/residual  value  *> 
begin 

k  :=  diecols  +  dj ; 

cpiz  HVectReadEl(pizel,k,rl) : 

if  regression  then 

kpiz  ;>  HVectReadEKresid.k.rl) 
else 

kpiz  HVectRead£l(pizel,k,rl) : 
end  {Replace  with  original  value} 
else  if  ((di  >  0)  and  (di  <  rovs-2))  and 
((dj  >  0)  and  (dj  <  col8-2))  then 
{*  Point  to  estimate  is  vitfain  borders  for  neighborhood  size  *} 
begin 

{«  Set  up  neighborhood  *} 
for  i  :«•  0  to  neighborhood- 1  do 
for  j  :=  0  to  neighborhood- 1  do 
begin 

kh  :=  (di-l+i)*col8  +  (dj-l+j); 
cn[i,j]  :»  HVectReadEl(pizel,kh,rl) ; 
if  regression  then 
kn[i,j]  HVectRead£l(re8id,kh,rl) 
else 

kn[i,j]  :=  HVectReadEl(pizel,kh,rl) ; 
end;  {for  j} 

{*  Estimate  pizel  value  «} 
if  (mi  <>  0)  and  (mj  <>  0)  then 
begin  {«  It’s  a  middle  position  e> 
cpiz  :=  (cn[0,0]+cn[0,3]+cn[3,0}+cn[3,3]) 

-  5*(cn[0,l]+cn[0,2]+cn[l,0]+cnCl,33+cn[2,0]  + 

cn  [2 , 3]  +cn  [3,1]  +cn  [3 , 2] 

-  5*(cn[l,l]+cn[l,2]+cn[2,l]+cn[2,2])) ; 
cpiz  :«  cpiz/64; 

kpiz  :=  0; 

for  i  :■  0  to  neighborhood- 1  do 
for  j  :»  0  to  neighborhood- 1  do 
kpiz  :=  kpiz  +  wMCi,  j]*kn[i,  j]  ; 
end  {*  Middle  *} 
else  if  (mj  <>  0)  then 

begin  {*  It’s  a  top  position  *} 
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cpix  :»  5*(cn[l,l]+cn[l,2])  -  cii[l,0]  -  cn[l,3]; 
cpix  :■  cpix/8; 
kpix  :=  0; 

for  i  0  to  neighborhood-l  do 
for  j  0  to  neighborhood-1  do 
kpix  :=  kpix  +  wTCi,  j]*kn[i,  j]  ; 
end  {*  Top  ♦} 
else 

begin  {•  It  nust  be  a  left  position 

cpix  :■  5*(cn[l,l]+cn[2,l])  -  cn[0,l]  -  cn[3,l]; 

cpix  cpix/8; 

kpix  0; 

for  i  0  to  neighborhood- 1  do 
for  j  0  to  neighborhood- 1  do 
kpix  kpix  +  wL[i, j]*kn[i, j] ; 
end;  {else}  (*  Left  *) 
end  {else  if} 

else  {♦♦♦  In  border  area  ♦**} 
begin 
cpix  :»  0; 
kpix  :=  0; 
end; 

{*  Add  regression  equation,  if  required  *} 
if  regression  then 
begin 

if  second_order  then 

reg  :=  c[0]  +  c[l]*sj/2  +  c[2]*si/2  +  c[3]*Sqr(sj/2)  +  c[4]*Sqr(8i/2) 
else 

reg  :■  c[0]  +  cCl]*8j/2  +  cC2]*8i/2; 
kpix  reg  -  kpix; 
end;  {if  regression} 

HNatHrtEl(CSOL,si,sj ,cpix,rl) ; 

HHatVrtEKKSOL.siiSj  ,kpix,rl) ; 
end;  {for  sj} 

HArrFree (pixel) ; 

HArrFree(resid) ; 

Urite( ’Calculating  statistics. . .  ’) ; 

Assign(fbo,fno+’ .OKI’) ; 

Reset (f bo) ; 

8un_c  :»  0;  stui_c2  :»  0;  bia8_c  :»  0; 

8un_k  :=  0;  sun_k2  0;  bias.k  0; 
n  0; 

for  si  0  to  rdin-l  do 
begin 
Mark_Ti*e ; 

for  sj  :«  0  to  cdin-1  do 
begin 

Read(fbo,pixel_o)  ; 

■i  si  nod  sub_sanple;  di  :•  si  div  sub.saaple; 

■j  sj  nod  sub.sanple;  dj  sj  div  sub.sanple; 

if  <<di  >  0)  and  (di  <  ro¥8-2))  and 

((dj  >  0)  and  (dj  <  col8-2))  and 

((■i  <>  0)  or  (nj  <>  0))  then 
begin 

pixel.c  :•  HMatReadEl(CS0L,si,sj,rl) ; 
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pixel.k  HMatReadEl(KS0L,ai,8j ,rl} ; 
n  :■  n  +  1; 

dit.c  (pixel_o  -  pixal.c); 
dil_k  (pixel_o  -  pixel.k) ; 
bias.c  :■  biaa.c  +  dif_c; 

bias.k  :■  bias_k  +  dif_k; 

siia.c  :■  sua_c  +  Abs(dif_c); 
sua_k  :■  sua.k  AbaCdif.k); 
sua_c2  :■  8ua_c2  +  Sqr(dif_c); 

8ua_k2  ;•  8ua_k2  *■  Sqr(dil_k); 

end;  {if} 
end;  {for  j} 
end;  {for  i} 

HArrFree(KSOL) ; 

HArrFree(CSOL) ; 
nean.c  ;■  sun_c/n; 

8tdT_c  :■  sqrt ( (n*aun_c2  -  Sqr(8iui_c))/(n*(n-l))) ; 

■ean_k  :»  aun_k/n; 

8tdv_k  :■»  8qrt ( (n*aiui_k2  -  Sqr(aum_k))/(n*(n-l))) ; 

z_value  :>■  (■ean_c-nean_k)/(8qrt((aqr(8td»_k)/n  +  8qr(8tdy_c)/n)))  ; 

Vriteln; 

Vriteln( ’Cubic  Error’); 

Writeln(’  Mean  •  ’ ,nean_c:5:2, ’  < ’ , (bia8_c/n) :5:2, ’) ’) ; 

VritelnC’  StDv  »  ’ ,8td»_c:5:2) ; 

Writeln(’Kriging  Error’); 

VritelnC’  Mean  ••  ’  ,8ean_k:5:2,  ’  ( ’ , (biaa.k/n) ;5:2, ’) ’)  ; 

VritelnC’  StDv  -  ’ ,8tdT_k:5:2) ; 

Vriteln; 

VritelnC ’Z  value  ■  ’ ,z_Talue:6:2) ; 

VritelnC’  n  ■  ’,n); 

EID. 


E.2  Quinn-Curtis  Routine  Description 

The  following  is  a  brief  desciption  of  the  Quinn  Curtis  procedures  used  in  the  program  above 

(22:). 

e  HMatDef/HVecDef  creates  a  matrix/vector  of  type  REAL. 

e  HMatWrtEl/HMatReadEl  writes  or  reads  one  element  from  a  real  matrix. 

e  HVecWrtEl/HMatReadEl  writes  or  reads  one  element  from  a  real  vector. 

•  HMultipleReg  fits  a  least  squares  multiple  regression  equation  to  a  data  set.  It  chooses  a 
model  which  minimizes  the  sum  of  squares  of  the  distances  between  the  observed  responses 
and  one  predicted  by  the  fitted  model. 
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•  HStatAnalysis  performs  a  statistical  analysis  of  regresssion  results  to  include  standard  er¬ 
ror,  multiple  correlation  coefficient,  coefficient  of  determination,  as  well  as  significance  of 
regression. 

•  HVectMean  determines  the  mean  of  a  vector. 

•  HVectStDev  determines  the  standard  deviation  of  a  vector. 

•  HArrFree  makes  available  memory  space  used  by  an  array. 

•  HGaussJordan  solves  a  system  of  linear  equations  using  Gauss-Jordan  method  with  partial 
pivoting.  Returns  the  inverse  of  the  orignial  matrix  as  well  as  its  determinant  and  solution. 
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